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Abstract 


Two methods for developing high order single step explicit algorithms on symmetric sten- 
<^s with data on only one time level are presented. Examples are given for the convection and 
lineariz^ Euler equations with up to the eighth order accuracy in both space and time in one 
space dimension, and up to the sixth in two space dimensions. The method of characteristics 
is generalized to nondiagonalizable hyperbolic systems by using exact local polynomial solu- 
tions of the system, and the resulting exEict propagator methods automatically incorporate 
the correct multidimensional wave propagation dynamics. Multivariate Taylor or Cauchy- 
Kowaleskaya eipansions are also used to develop algorithms. Both of these methods can be 
applied to obtain algorithms of arbitrarily high order for hyperbolic systems in multiple space 
dimerous. Cross derivatives are included in the local approximations used to develop the 
algorithms in this paper in order to obtain high order accuracy, and improved isoptropy and 
stability. Efficiency in meeting global error boimds is an important criterion for evaluating 
algorithms, and the higher order algorithms are shown to be up to several orders of magnitude 
more efficient even though they are more complex. Stable high order boundary conditions for 
the linearized Euler equations are developed in one space dimension, and demonstrated in two 
space dimensions. 



I: Introduction 


Hyperbolic systems include familiar examples such as the linear systems for convective 
tramsport, acoustics, and electromagnetics, and the nonlinear Euler equations of fluid mechan- 
ics ([18], [39]). These examples have very practical applications in areas such as aircraft noise 
([24], [33], [35]). Numerical methods for hyperbolic systems have a broad history ([8], [20], [36]). 
There is an increasing interest in the use of computational fluid dynamics techniques for these 
systems ([26], [27], [32], [34], [42]), but with imusually severe accuracy requirements ([16], [25]). 
Algorithms are needed for hyperbolic systems that are able to propagate a wide remge of 
wavelengths with high accuracy over long distances. 

Compact difference methods ([5], [29], [37]) are generally viewed as being highly accurate, 
and can be tailored to produce both high order accuracy and high resolution [22], which is 
defined as the ability to propagate relatively high frequency waves with the correct velocity. 
Compact difference methods generally use separate treatments for space and time, either with 
different [40] or the same [14] order of accuracy. The dispersion relation preserving method 
[38] is similar to compact difference methods in its spatied treatment, but it addresses the 
relationship between space and time treatments by using some of the degrees of freedom of 
temporal data that it requires in order to reduce the dispersion or phase speed errors in 
each time step. The dispersion relation preserving scheme provides fourth order accuracy in 
space emd second or third order in time on a seven point stencil with data from four time 
levels. High resolution is frequently stated in terms of the number of grid points required 
to accurately propagate a normal mode, or algorithm performance relative to a grid scale, 
which does not address the question of efficiency with respect to meeting a stated global error 
boimd. Compact difference methods do not provide a genereil efficient high order time stepping 
method. 

A variety of high order finite difference methods use multiple time steps. The dissipative 
two-four method [13] is a two step generalization of the Lax-Wendroff method, with second 
order accuracy in time and fourth order accuracy in space. The two-four method has versions 
([1],[13]) that are similsu' to the MacCormack [28] operator splitting method. The third order 
difference method of Bmstein, Mirin and Rusanov ([3], [31]) uses intermediate time steps and 
a correction to obtain third order accuracy in space and time on a five point stencil. The 
modified equation approach [6] uses data on three time levels to obtain fourth order accuracy 
in space and time. A five-six finite difference scheme has been developed [43] on a seven 
point stencil in one dimension with a six stage time marching method, providing fifth order 
accuracy in time and sixth order in space. This method has an optimized variation which 
relaxes its order of accuracy in order to increase its resolution. The particular approach of 
directional sphtting in multiple space dimensions violates the physics of nondiagonalizable 
hyperbolic sytems in a fimdamental way, and in generad, the use of multiple time steps raises 
issues such as intermediate time level boundary conditions, £trtifical dissipation, significant 
time step constraints, starting values, and efficiency. 

Many classical approaches to developing numerical methods can be viewed as simtdtane- 
ously treating the spatial representation of data and the temporal evolution of the system. Lax 
Wendroff [21] methods can be viewed as using a second order Taylor series expansion in time. 
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with time and space derivatives related by the partial differential equation. Semi Lagrangi^ 
methods use the method of characteristics [41] and incorporate the geometric behavior of the 
solution in space and time. Godunov [10] methods use the solution of a Riemann problem to 
approximate the physics of shocks within a control volume in spax:e and time. Fimte volume 
methods use integral forms of conservation laws in space and time, and are close in spint to the 
general method of control volume analysis. The finite analytical method [4] also combines the 
treatments of space and time by using a local Fourier decomposition m space and a sep^ation 
of variables analytical treatment for time evolution. These approaches generally are limi e 

in either accuracy or applicability. , . , i i * 

This paper presents and compares two different approaches to algorithm developmen 
for hnear hyperbolic systems in multiple space dimensions, the use of local exact polynomial 
solutions in space and time, or exact propagators, and the use of multivariate T^lor senes, 
or Cauchy-Kowaleskaya expansions. Local exact polynomial solutions can be obtained tor 
diagonahzable hyperbolic systems by the method of characteristics, and for nondiagonaliz- 
able systems in multiple space dimensions by requiring that a general polynom^ expansion 
in space and time exactly solves the partial differential equations, with all of the unknown 
expansion coefficients that involve time expressed in terms of spatial coefficients. The Oauchy- 
Kowaleskaya procedure ([9],[17]) is an equivalent method for obtaming time expansion coef- 
ficients. Local exact solutions which incorporate the multidimensional wave dynamics of the 
hyperbohc system can be viewed as a correct way to extend the method of charactenstics to 
nondiagonalizable systems. If the Taylor series methods are considered only at the spatial 
center of the local expansions, then they can be viewed as a Taylor senes expansion m time 
alone, similar to the second order Lax Wendroff method [21]. The exact propagator and mul- 
tivariate Taylor series approaches produce single step exphcit algonthms that have the same 
order of accuracy in both space and time, while using data from only one time level, and that 
can be extended to arbitrarily high order and multiple space dimensions for nondiagonalizable 
systems. These two approaches to algorithm development produce the same methods in one 
space dimension, and different methods in multiple dimensions. 

The exact propagator and multivariate Taylor series approaches to algorithm development 
can both be viewed as direct applications of the general use of multivariate approximating 
polynomials in space and time ([2], [8]). As a related example of this te^nique, in a dis- 
cussion from the Taylor series perspective about extending Leith’s method from one o wo 
space dimensions, Roache [30] makes the point that a scheme in two space dimensions would 
have to include cross derivative terms in order to obtain high order accuracy. The use of cross 
derivative spatial terms has been shown to also improve isotropy and stabihty, for example in 
the development of finite volume methods for multidimensional advection [23], and in the de- 
velopment of a nearly exact second order algorithm on three time levels for the wave equation 
[7]. Including sufficient cross derivative terms in a multidimensional algorithm is required for 
high order accuracy, and improved isotropy and stability, but it does not necessarily produce 
a locally exact solution with the correct multidimensional wave dynamics. A fundamental 
viewpoint in this paper is to approximate the solution of a system of partial <fifferential equa- 
tions as a whole, instead of approximating separate derivative terms in particular equations. 
This viewpoint is expressed in the two general approaches to algorithm devebpment that are 
used throughout this paper, which in turn are realized in a sequence of particular algorithms 
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for linear hyperbolic systems in one and two space dimensions. 

Development of high order algorithms for the two dimensional linearized Euler equations is 
a prominant goal of this paper, but various issues are addressed in the simplest possible context 
by developing algorithm for other systems. Examples for the convection and linearized Euler 
equations will be shown in one space dimension up to eighth order in both space and time, and 
in two spzwre dimensions up to sixth order in both space and time. In the second section of this 
paper, a fomi;h order algorithm is presented for the scalar first order linear wave equation in 
one space dimension, and is used to relate several viewpoints of algorithm development. In the 
third section, algorithms for the one dimensional Hneanzed Euler equations are presented, with 
second, fourth, sixth and eighth order accuracy in space and time. The accuracy and relative 
efficiency of these algorithms is examined for short and long time calculations, with high and 
low error boimds [19]. The problem of creating stable high order boundary algorithms for this 
system is also treated. In the fourth section, algorithms are developed on symmetric stencils 
for the convection equation in two space dimensions with second, fourth and sixth order 
accuracy in space and time. Relative efficiency, stencil choice, and the differences between 
exact propagator and Taylor expansion algorithms are examined. In the fifth section, exact 
propagator and Taylor series algorithms with second, foiurth, and sixth order accurcy in both 
space and time are developed for the linearized Euler equations in two space dimensions. A 
method for developing local exact polynomial solutions to the two dimensional linearized Euler 
equations is described and used for developing exact propagator algorithms. This method for 
developing local exact polynomial solutions is a departure from the method of characteristics 
which is used throughout the rest of the paper. The relative effect of algorithm type and 
order is compared, and a high order implementation of a new imobtrusive outflow boundary 
condition by Hagstrom [15] is demonstrated. The sixth section is a brief concluding summary. 
Several longer formulas have been included in appendices. 
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II: A Fourth Order Algorithm For The ID Convection Equation 


Consider the scalar first order linear wave equation in one space dimension for u{x,t), 


du ..du 


( 1 ) 


where M is the constant mean convection speed. The initial value problem for this equation 
with u(i,0) = «i(x) for I € can be immediately solved by the method of characteristics, 
and u(x,t) = ui{x-Mt) for x G and 0 < t. With this weU understood ex^ple it is possible 
to quickly develop high order edgorithms, and to readily present an algorithm from various 

A fourth order numericed algorithm can be developed for this problem by the method of 
charateristics with a local quartic interpolation to u at time on a five point central stencil. 
A uniform grid is used, with (xi,t„) = (z7i,nfc) for integer i and n, where h = Ax and k = At 
are the uniform mesh spacings in x and t. The quartic spatial approximation to u around x^ 


at tn can be written in local coordinates as 


u(Xi + X^tn) ^ Ua(x) = Uo + + U2X^ + U3X + U4X , 


where the coefficients are not indexed with respect to the mesh point. The method of un- 
determined coefficients and the known data on the five point stencil immediately yield the 
xinknown coefficient solutions 


Uo = Ui 


Ui 


y>2 




U4 


= T^(“r-2 - 8u?_. + 8u”+, - «r+2), 


I2h 




2Ah^ 


( 2 ) 


= + 2“?-i - 2«?+i + <‘"+ 2 ), 


12/l3 


= “”+2)- 


2Ah^ 


Notice that ~ spatial approximation to u about Xi at t„ is 

essentially a truncated Taylor series in x. The method of characteristics and the local solution 
approximation at time tn now give 


u{Xi,tn+l) = u{Xi - Mk,tn) « Ua(~Mk) = Y, = < 


n+1 


a=0 


where is implicitly defined. This algorithm correctly incorporates the wave dynamics of 

equation (1) by using the method of charateristics. Notice that the time propagation is exact, 
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so that the error at each time step is due to the local spatial interpolation ua. This algorithm 
can be rewritten in the familiar form of a conventional single step explicit finite difference 
method, 


u 


n+l 


= E 


i-fr? 


r=-2 


where if the CFL number is A = then 


a+2 = 

—(+2 - A - 2A^ 

+ A^), 

0+1 = 

A(_4 + 4A + a2- 
0 

-A^), 

oq = 

i(4 - 5A^ + A*). 


a_i = 

j(+4 + 4A - A^ - 

D 

-A^ 

a_2 = 

A(-2 - a + 2A= 

+ A3). 


In this form the algorithm requires five multiplies and four adds at each grid point for each 
time step. Standard analysis shows that this method is fourth order accurate in both space 
and time, with truncation error 


-{u(Xi,tn+l) - 


)='-‘(-j^)(1-A^)(4-A^) 




+0[h% 




d^u 

dx^ 


The leading order error term is dispersive, which is natural for algorithms with symmetric 
central stencils. Numerical experiments suggest that the method is stable for |A| < 1. 

If the local expansion coefiicients axe interpreted in terms of space derivatives, and if u is 
an exzw:t solution to equation (1), then 


u 


n+l 


= E 


Q=0 


{-Mk)° d^ua 
o! dx° 




d°‘u 


cr=0 


a! 


dx 




a=0 


a! dt 


d^u. 

a V^t 9 


This form of the algorithm is a fourth order expansion in the time step size ifc, and is similar to 
the standard second order Lax-WendrofF method [21] expressed as a tnmcated Taylor series in 
time. If the general form of the solution to equation (1) is used, then u can be approximated 
in both local coordinates x and t near (x,,<„), with u(xi + x,<„ + 1) sa ua{x - Mt). If the 
interpretation of the expansion coefficients axe used, then this local approximation to « in x 
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and t can be written as 


u(Xj + x,tn+ t) «ua(x - Mt) 

4 

= - Mty 


a=0 

4 




A 

a=0 
4 a 


a=0 B=0 ^ ' 


where = ua{-Mk). This form of the solution is actually just a truncated Taylor series 

in both space and time simultaneously. This truncated Taylor polynomial is also an exact 
polynomial solution to equation (1). This particular exact local solution approximates the 
solution u for points {x,t) with a domain of dependence contained in the interval [xj_ 2 ,x,+ 2 ]- 
There are four related interpretations of this fourth order algorithm: a locally exact solu- 
tion derived from the geometric method of characteristics with a local truncated Taylor series 
approximation in space for its initial data; a conventional finite difference method; a truncated 
Taylor series expansion in time; and a locally exact polynomial solution derived analytically 
from a truncated Cauchy-Kowaleskaya or bivariate Taylor series expansion in space and time. 
The first interpretation ensures that the numerical method properly represents the wave dy- 
namics of the partial differential equation, since the characteristic behavior of the equation is 
incorporated into the numerical solution. The second and third interpretations ground this 
algorithm in the mature tradition of finite difference methods. The fourth interpretation pro- 
vides an avenue for generalization and extension, since a locally valid exact analytic solution 
is possible for hyperbolic systems in higher space dimensions. Notice that no consideration 
has been directly given to the problem of finite difference approximation to any derivative, 
but rather to the interpolation of the known data on a stencil, and to an exact local solution 
to the partial differential equation. A key shift in perspective is away from the detmls of 
approximating separate terms in an equation, and toward the approximation of a solution to 
the equation. 
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Ill: Numerical Algorithms For ID Linearized Euler Equations 


Consider now the linearized Euler equations in one space dimension, in particular the 
isentropic case in the form of the nondimensionalized system 

du du dv 

dp dp du ^ ^ 

where M is the constant mean convection speed in terms of Mach number, and where p and u 
are the pressure and velocity of the distiurbance. A general discussion of the linearized Euler 
equations can be found in Kreiss and Lorenz [18]. System (3) can be diagonalized and written 
with Riemann variables in the equivalent decoupled form 

d<jJt , . dujA 

M-1 ^=0, 

dt dx 

where = |(u — p) and u ?2 = ^(u + p). The initial value problem for u and p with u(x,0) = 
«*(x) and p(x, 0) = pt(x) for x 6 3? is equivalent to a corresponding problem for ui and 
u }2 with a;i(x,0) = |(ui(x) — pi(x)) and u>2(a;,0) = |(ui(x) + pi(x)). Each equation in the 
Riemann variable system (4) can be solved separately by the method of characteristics, md 
these solutions can be used to obtain the general solution for system (3), 

«(x, t) =^(ui(x — (M + l)t) + pi(x — (M + l)t)) 

+i(iii(x — (M — l)t) — pz'(x — (M — 1)<)), 

1 ( 5 ) 

=2(“*(® ~ (^ + l)t)+pi(x — (M + 1)<)) 

— ^(ui(x — (M — 1)<) — pz(x — (M — l)t)). 

Note that this system has two characteristics with separate imidirectionzil solutions that are 
constant along their characteristics, and that are travelling with the characteristic velocities 
M — 1 and M + 1, respectively. 

III-A: A General Algorithm Development 


A general development can simultaneously provide algorithms of the second, fourth, sixth 
and eighth orders. A local polynomial interpolation to u and p in x about (x,-,fn) can be 
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written in the general form 


Or 

u(Xi + X,tn) « ua(x) — ^ UaX", 

a=0 

Or 

p(xi + X, t„) « pa(x) = ^ Pqx“, 

0=0 


( 6 ) 


where indexing with respect to the mesh point is suppressed, and where Or = 2, 4, 6 , or 8 is 
the order of the interpolation, on central stencils of from three to nine points. The method 
of tmdetermined coefficients using the data on the stencils immediately yields the unknown 
coefficients in terms of the known data. The coefficients have standard central finite difference 
forms, such as for the quartic case in (2). The expansion coefficients can be interpreted in 
terms of spatial derivatives, with 


Ua = 
Pa = 


1 d°ua 

1 d°pa 
a! dx°‘ 


( 0 )~ 


( 0 ) « 


1 

a! dx°‘ 
1 5“p 
a! dx^ 


(x j, ), 


5 tn). 


A new solution value at (xj, ) is obtained from the general solution form (5) with the local 
spatial interpolation ( 6 ) as initial data. 


u(xi, =^(ua(-(M + 1)A:) + pa(-(M + l)k)) 
+ hua(-(M - l)Ar) - pa(-(M - l)fc)), 
p{xi,tn+l) ~ p”'*’^ + 1)^) +pa{ — (M + l)fc)) 

- l)k) - pa(-(M - l)k)). 


Algorithm form (7) is obtained by applying an exact propagator to a local polynomial inter- 
polant, so that form (7) correctly incorporates the mtdtiple characteristic wave dynamics of 
equation (3). Algorithm (7) can also be arranged as a truncated time series expansion in k. 
On the nine point central stencil the eilgorithm for u can be rewritten as 

= uo — k(pi + Mui) + fc^( 2 Mp 2 + (M^ + 1 )« 2 ) 

- P((l + 3M^)p3 + M(3 + M^)u3) 

+ fc^(4M(l + M^)p4 + (1 + 6 M^ + M*)u4) 

+ jfc®(-(l + lOM^ + 5M^)ps - M(5 + lOM^ + M‘‘)us) ( 8 a) 

+ fc®(M (6 + 20M^ + 6M*)pe + (1 + 

+ 4 . 2iAf2 + 35M^ + 7 M®)p7 - M(7 + 35M^ + 21M^ + M^)ur) 

+ ib®(M (8 + 56M^ + 56M^ + 8 M®)ps + (1 + 28M^ + 70M* + 28M® + M^)us), 
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and for p as 

=Po- Hui + Mpi ) + P(2Mu 2 + (M^ + l)p2) 

- + 3M2)u3 + M(3 + M^)p3) 

+ A:^(4M(1 + M2)u4 + (1 + 6M2 + M*)p^) 

+ A:®(-(1 + lOM^ + 5M^)us - M(5 + lOM^ + M“)p5) (86) 

+ k^{M{6 + 20M^ + m*)u6 + (1 + ISM^ + 15M^ + M®)pe) 

+ *^(-(1 + 21M^ + 35M‘‘ + 7M®)«7 - M(7 + 35M^ + 21M^ + M^)pj) 

+ fc®(M(8 + 56M^ + 56M^ + 8M®)us + (1+ 28M^ + 70M^ + 28M® + M®)p8)- 

Algorithm form (8) can be used to obtain the time series expansion forms of the second, fourth 
and sixth order algorithms, by simply taking the terms up to the appropriate order in k. Note 
that the symmetry of (3) in u and p is reflected in (8). The familiar grid ratio A = |^ is implicit 
in algorithm form (8), since each expansion coefficient Uq or pa is a finite difference form with 
in its denominator. A general finite difference form can be derived from (8), with 

R 

= E 

R 

= E 

r=-R 

where the coeflBcients axe polynomials in A and M, and where iZ = 1, 2, 3, or 4 depending upon 
which stencil is being used. The coefficients and br are used symmetrically with respect to 
u and p in algorithm form (9), and reflect the symmetry of the equation (3) in u and p. 

Algorithm (7-9) on a three point central stencil has Or = 2 in (6) and i? = 1 in (9), 
and represents just the first line of equations (8a) and (8b). The truncation error for this 
algorithm is 


1/ / . \ n4-l^ ^ -.^du dp 


+ - (3 + - (1 + 


d^u 


dx 


d^p 

dx^ 


+ 0[h% 




at ax ox 


+ )(1 - (3 + + A"(g)(l - (1 + 

+ 0[h% 


Clearly, this method is consistent with equations (3), it is second order accurate in space and 
time, and it is dispersive. Algorithm (7-9) on a five point centrzd stencil has Or = 4 in (6) 
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lines of equations (8a) and (8b). This 


and R = 2 in (9), and represents just the first three 
algorithm has the tnmcation error 


1 , , du » , du dp 

-(w(xi,<„+i)-< ) = ^ 

+ h\=^ + M{3 + M^)^- M(5 + lOM^ + 

+ + (1 + 3M^)^ - (1 + lOAf^ + 

+ OlA'l, 

1 dp dp du 

+ '■‘(^ + M(3 + - M(S + 


30 

+ A <(^ + (1 + - (1 + lOM ^ + 8 M ‘)^)0 

+ 0[h% 


This fourth order accurate method is clearly consistent with (3) and dispersive. The truncation 
errors for the seven and nine point central stencil methods also verify that these methods are 

sixth and eighth order accurate, respectively. ^ 

The algorithms in this section can be viewed as characteristic or semi Lagrangian meth- 
ods with a local polynomial spatial interpolation, as truncated Taylor series in time, or as 
conventional explicit finite difference methods. The exact solution form and the mterpreta- 
tion of the local polynomial expansion coefficients can be used to provide a local polynonual 
solution to the differential equation that can be interpreted as a midtivariate Taylor series or 
Cauchy-Kowaleskaya expansion in space and time. The algorithms in this section are derived 
from an exact solution to the partial differential equation, and find a new solution value by 
using the correct wave dynamics of a local spatial interpolant. 


III-B: Numerical Comparisons 

Numerical tests of these four algorithms for equation (3) are obtained from an imti^ 
value problem with data u(x,0) = 0 and p(x,0) = Sin(Trx), for —1 < x ^ 

boundaries at x = ±1. The exact solution for this problem can be obtained from (5). Table 
A gives data for this problem with the four algorithms atM = 0andA = ^ = 0.8, on a senes 
of mesh sizes, and calculating out to the fixed time t = 10, which corresponds to five periods 
of wave propagation for the given initial data. In Table A, ^ is the n^ber of grid points 
in [—1, 1], or per wavelength for the initial data, nio is the number of time steps required to 
reach t = 10 with A = 0.8 at the given grid resolution, and the columns are identified by the 
order Or of the algorithm which produced the data in the column. The error data in Table A 
is the TTiavlTrmTn absolute accumulated error in u or p for x e [-1, 1] at t = 10, so it reflects 

10 


I 



the error of the algorithms in both space and time. The mesh sizes that eure used are for mesh 
refinements by successive halving of the grid size. 


Table A: Data From The Four Algorithms For The ID LEE 



u(x,0) 

= 0 and p(x, 0) = 

= Sin(7ro:), A 

= 0.8,M = 

:0 



Maximum Error 

in « or p at < = 10 


2 

h 

nio 

Or = 2 

Or = 4 

Or = 6 

Or = 8 

4 

25 

1.035D+0 

7.885D-01 

4.253D-01 

2.078D-01 

8 

50 

6.760D-01 

9.141D-02 

l.llOD-02 

1.387D-03 

16 

100 

2.586D-01 

7.122D-03 

2.158D-04 

7.000D-06 

32 

200 

7.133D-02 

4.644D-04 

3.551D-06 

2.910D-08 

64 

400 

1.811D-02 

2.932D-05 

5.620D-08 

1.157D-10 

128 

800 

4.539D-03 

1.837D-06 

8.808D-10 

3.438D-13 

256 

1600 

1.135D-03 

1.149D-07 

1.283D-11 

9.484D-13 

512 

3200 

2.839D-04 

7.182D-09 

6.303D-13 


1024 

6400 

7.097D-05 

4.527D-10 

3.776D-12 



The orders of accuracy of the algorithms can be verified by calculating the factor by which 
the error decreases with the grid resolution. The data from the method with Or = 2 at grid 
sizes ^ = 512 and ^ = 1024 gives, 


1 2.839 X lO-'* 

Log[2r^^h.097 x IQ-^ 


] = 2 . 00 , 


so that the error decreases by a factor of 2^ “° when the mesh size is halved, and the method 
as implemented actually is second order in both space and time. For the method with Or = 6, 
the error decrease from = 128 to ^ = 256 is by a factor of 2®*^®, but the decrease from 
j = 256 to = 512 is only by a factor of 2*-^®, while the error actually increases from ^ = 512 
to ^ = 1024. The differences between these estimates of the order of accuracy for this method 
is due to the contamination of the finer grid results by roimdoff error at the resolution limits 
for these calculations in double precision. The data from j = 128 and ^ = 256 can be 
accepted as showing sixth order accuracy for this method. A similar effect can be seen in the 
data for the method with Or = 8, but the order of the method can be obtained from the data 
at = 64 zuid ^ = 128, which shows an error decrease by a factor of 2*’^® and eighth order 
accuracy. The order of accuracy of a finite difference method is an asymptotic concept that 
apphes in the limit as the mesh size converges to zero, but due to roundoff errors, niunerical 
computations in a given precision to a fixed simulation time with a particular algorithm have 
an effective lower boimd on usable grid resolution that is finite, even though the bound can 
be affected by programming practices. 

The second order algorithm with \ ^ has a maximiun absolute error in u and p at 

t = 10 of 1.03451. At t = 10, the numerical solution from this simulation shows a maxdmmn 
absolute value of 1.55885 x 10"^ for u, and 3.45145 x 10“^ for p, while the exact solution 
has u identically 0, and extreme values of ±1 for p. The larger than one absolute error arises 
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because the p solution is a half wavelength out of phase, which is possible even at a short time 
with a very coarse resolution and a dispersive method. Notice in Table A that the maadmum 
errors at f = 4 range from 0[1] to O[0.1], decreasing slightly as the order of accuracy of the 
method increases. This shows no outstanding advantage for the higher order methods at this 
coarse grid resolution. The actual advantages of higher order methods only become apparent 
when the qualities of the methods are not confounded with marginal data resolution. 

Table B presents data from the same periodic problem as in Table A, but at f = 1000. In 
Table B the number of time steps is now given as riiooo, the number required to compute to 
t = 1000 with A = 0.8 and the given grid resolution. The errors in Table B from the second 
order method at | = 16 and the fourth order method at f = 8 are produced by a very large 
phase error, just as the error in Table A for the second order method at ^ = 4. In general, 
the errors in Table B tend to be larger than the corresponding errors in Table A by a factor 
of 100, which is the ratio of niooo to nio at any given grid resolution. These ^gorithms are 
derived with an exact propagator for equation (3), so that the discrete approximation of the 
time evolution can be expected to have 0[1] eigenvalues and a linear growth trend in error. 
But notice that the maximum absolute errors on the coarsest grids never become much larger 
than one, or that the error growth has an asymptotic limit as the simulation time increases, 
since the errors eventually become of the order of the solution. Also notice that the noise level 
where accumulated roundoff error begins to significantly affect the solution has increased from 
O[10“^^] at f = 10 to O[10"®] at t = 1000. 


Table B: Data From The Four Algorithms For The ID LEE 
u(i,0) = 0 and p(a:,0) = Sin(Trx), A = 0.8, M = 0 
Maximum Error in tt or p at i = 1000 


2 

h 

4 

8 

16 

32 

64 

128 

256 

512 

1024 


UlOOO 

2500 

5000 

10000 

20000 

40000 

80000 

160000 

320000 

640000 


Or = 2 Or = 4 Or = 6 Or = 8 


l.OOOD+0 

l.OOOD+0 

l.OOlD+0 

7.401D-01 

1.215D+0 

4.327D-01 

1.131D-01 

2.838D-02 

7.096D-03 


l.OOOD+0 

1.004D+0 

5.407D-01 

4.614D-02 

2.931D-03 

1.837D-04 

1.150D-05 

7.000D-07 

1.294D-08 


l.OOOD+0 

7.787D-01 

2.145D-02 

3.552D-04 

5.618D-06 

9.263D-08 

9.364D-09 

1.816D-08 


l.OOOD+0 

1.265D-01 

6.998D-04 

2.909D-06 

9.548D-09 

4.590D-09 

7.987D-09 


If interpolation is used with the data in Table A, then this data suggests that in order 
to achieve a maximum absolute error of less than 5.0 x 10 at t = 10, the second order 
method requires ^ > 215, the fourth requires ^ > 16, the sixth ^ ^ 8, and the eighth ^ > 7. 
For a grid ratio of A = the number of time steps nio required to compute to t = 10 is 

= IM. The total number of multiplications required for each computation on the domain 
— 1 < X < 1, with j mesh points per unit interval, using a stencil with n, grid points, with 
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two equations using data from two functions, and for nio time steps is 


toted mtiltiplications = nio x 2h x n, x 4 = 


SOh^rig 

A 


With modem RISC processors, a single cycle multiply emd eidd instruction is possible, so that 
additions can be neglected in floating point operation counts. In order to meet the 5 x 10“^ 
error bound at the time t = 10, the ratio between the number of multiplications required by 
the second and fomrth order methods is 


2152 X 3 

162 X 5 


108.34 


it is 2.86 for the fourth to the sixth order, and 1.02 for the sixth to the eighth order. If the 
error boimd is decreased to 5 x 10~® at < = 10, then the second order method requires ^ > 673, 
the fourth ^ > 31, the sixth ^ > 14, and the eighth ^ > 8. In this case, the ratio of the 
nmnber of multiplications required by the different methods is approximately 282.79 for the 
second to the fourth order, 3.50 for the fourth to the sixth order, and 2.38 for the sixth to the 
eighth order. If the error bound is kept at 5 x 10~^ but the time is increased to < = 1000, then 
the second order method requires j > 1989, the fourth ^ > 60, the sixth ^ > 16, and the 
eighth ^ > 10. In this case, the ratio of the nmnber of multiplications required by the different 
methods is approximately 659.4 for the second to the fourth order, 10.04 for the fourth to the 
sixth order, and 1.99 for the sixth to the eighth order. 

The difference between low and high order algorithms is not very significant if the data 
is not resolved well enough for any of them to be able to accmately propagate meaningful 
information to a useful time. It is significant that small errors can be obtained at relatively long 
times, even though there is a limit to what can be achieved with double precision arithmetic. 
For example, if the periodic problem is computed with the eighth order method at ^ = 16 
out to t = 100000, or 50000 periods, then the maximum absolute error 3.022D-04 is invisible. 
The higher order methods axe clearly more eflicient in terms of the total number of floating 
point operations required to compute to a fixed simulation time with a given error bound, 
and the comparative efficiency of the higher order methods increases as the simulation time 
is increased or as the error bound is decreased. There is a particularly great advantage shown 
by the fourth order method compared to the second. These observations should be expected, 
since they are predicted by the analysis of Kreiss and Oliger [19]. Higher order methods can 
be more complex, requiring additional research and development time as well as extra care 
in implementation, but the value of the efficiency of higher order methods increases either if 
they are incorporated in codes that are frequently used, or if specific computations are made 
possible only with the use of a highly efficient algorithm. 

III-C: High Order Boundary Conditions 


Propagation algorithms must be complemented by high order boimdary algorithms in 
order to be useful. Stable high order boimdary algorithms can be developed by using the 
exact propagator approach. The general algorithm development for equation (3) uses a local 
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spatial interpolant (6) at time t„ as initial data, and the general solution form (7), to obtain 
an exact solution to the linearized Euler equations which correctly incorporates the wave 
dynamics of the local interpolant. This exact solution can be written in local coordinates as 

u(xi + + i) » ua{x,t) 

= kj] «a(x - (M + l)f)“ + Pa{x - (M + 

^ a =0 cr =0 


+ l(f] U„(l - (M - !)()“ - ^ Va(T -(M- 1)0°). 

^ a =0 cy =0 

p(xi + x,t„ + t) « pa{x,t) 

Or Or 

= Ua{x - (M + l)t)“ + Y, Pa{x - (M + l)t)“) 

^ a =0 «=0 


( 10 ) 


- - (JW - ')*)“ - -(M- !)()“). 

0^=0 of=0 

where Or is the order of the interpolant, and where and Pa are interpolation coefficients 
obtained by the method of undetermined coefficients. Note in (10) that an extra free variable 
has been added to the domain of ua and pa from (6), and that the general numerical algorithm 
in either form (7) or (8) can be obtained from (10) as = ua(0, k) and p”"*” = pa(0, k). 

The local exact solution (10) can be interpreted as a multivariate Taylor series or Cauchy- 
Kowaleskaya expansion in space and time. The local exact solution (10) represents more than 
just an approach to obtaining a numerical approximation to u and p at (xj,t„+i), it is also 
a valid representation of the evolution in space and time of the local spatial interpolant (6) 
at all points which have their domain of dependence entirely within the interpolation stencil. 
The Riemann variables uJi = ^(u — p) and u >2 = ^(w + p) can also be locally approximated in 
space and time, with woi = |(ua - pa) and ua 2 = \{ua + pa). For subsonic mean flow with 
going Riemann variable u}\ will have a valid representation up to the left end 
of the stencil, and the right going Riemann variable 013 will have a valid representation up to 
the right end of the stencil. The general boundary algorithm idea is that at the boundaries of 
the computational domain, the outgoing Riemann variables are calculated, and appropriate 
boimdary conditions are used to provide enough degrees of freedom of information to determine 
the primitive variable solutions at the boimdary. If there are any points between the center and 
boundary points of a boundary stencil, then solution values are computed at these intermediate 
points with (10). All of these different calculations use a single spatial interpolation on one 
boundary stencil, and a consistent representation of the local evolution of the interpolant over 
the interval from the stencil center to the boimdary of the computational domain. 

In the case of Or = 2 on a three point stencil, in addition to the new values at the stencil 
center the only calculation that is needed is the outgoing Riemann variable. At the stencil 
center, either (8) or (10) give 

= uo — k{pi + Muj) + P(2Mp2 + (M^ + l)«2)j ( 11 a) 

p""*"^ = po — k{u\ + Mp\) + k^{2Mu2 + (Af^ + 1)P2)* 
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The left going Riemann variable is obtained from (10) at as 


=2((“o -po)-h(«i -Pi) + h^(«2 -P2)) 
^ -pi) - 2 h(u 2 -P2)) 


2 

(1 - Mf 


(«2 -P2), 


( 116 ) 


and the right going Riemann variable at (xj+i,<„+i) as 


‘*^2”+!^ — Po ) + ^(“1 + Pi) + h^(u2 +P2)) 

-—^—((“1 +Pi) + 2 h («2 +P2)) 


(1 + My 


(U2 +P2), 


(11c) 


where the interpolation coefficients for u are 


uo = ui = - <- 1 ), and «2 = §^(“”+1 “ 2«r + “”- 1 ), 

with similar forms for p. H the three point stencil is at the left of the computational domain, 
then (11a) and (11b) are used, and if the three point stencil is at the right of the computational 
domain, then (Ha) and (11c) are used. In either case, a single interpolation stencil is used 
to provide the coefficients and Pa for all three calctdations, and an additional degree of 
freedom of information is required to obtain the primitive variable solutions at the boundziry 
point. 

In the case of Or = 4 on a five point stencil, the required calculations are for u and p 
values at the stencil center and one intermediate point, and for one outgoing Riemann variable 
at a stencil edge. At the stencil center, either (8) or (10) give 


=uo - k(pi + Mui) + k^(2Mp2 + (M^ + l)u2) 

- k^((l + 3 M^)p3 + M(3 + M^)u3) 

+ k*{AM{\ + M^)p4 + (1 + 6M^ + M^)u4), 
p|*+i = po - k{m + Mpi) + k'^{2Mu2 + (Af2 + l)p2) 

- ^"*((1 + 3A/2)u3 + M(3 + M2)p3) 

+ jfc^(4M(l + M2)«4 + (1 + 6M^ + M'*)p4), 


where the coefficients for u are given by (2), with similar forms for p. In a boundary stencil at 
the left of the computationed domain, the intermediate solution values between the cell center 
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zind the boiindaxy point are given by (10) as 

=Uq — hu\ + h^U2 — + h^Ui 

-\-k[M{—ui + 2hu2 — Zh?u^ + 4/i^U4) 

+ {—pi + 2/ip2 — Zh^p3 + Ah^p^)) 

+A;2((1 + M^)(u2 - Zhuz + 6h^Ui) 

+ 2M{p 2 - Zhp3 + Gh'^P4)) 

+k^{M{3 + M^)(-U3 + 4hu4) 

+ (1+ 3M^)(-P3 + 4hp4)), 

p^^i =Po — hpi + h^P 2 — h^Ps + h*p4 

+k{M{-pi + 2hp2 - Zh'^pz + 4h^P4) 

-j- ( — til "i" 2hu2 — Zh^U3 + 4/1^114)) 

+ k\{l + M^)(P2 - Zhp3 + 6h^P4) 

+ 2M{u2 — Zhuz + 6h^U4)) 

+k\M{Z + M^){-P3 + 4hp4) 

+ (1 + 3il/^)(— U3 + 4/1U4)), 

and the botindary point value for the outgoing Riemann variable is given by (10) as 

u?i|*+2^ =^Wi,o - ^1,1 + - 4h^^i,3 + 

H“fc(l — Af)(— — 2hoji^2 + — 16h^coi^^) 

+fc^(l — M)^(^wi ,2 — 3/iu;i,3 + 12 /i^u;i, 4) 

+k^{l - - 4/iu;i,4) 

+fc"(l-M)"^a;i,4, 

where ~ 5(^0 — Pq')i the coefficients Uq, and pa given by (2). 

In the case of Or = 6 on a seven point stencil, u and p values are needed at the stencil 
center and two intermediate points, and one outgoing Riemann variable is needed at a stencil 
edge. If the three botmdary treatments for Or = 2, 4 or 6 are used with Dirichlet bounday 
data for u and p, then they maintain the order of accuracy of the propagation algorithm m 
both space and time, and they are stable if \M^\ < 1. In the case of Or = 8 on a nine 
point stencil, if this general procedure is followed, with u and p values calculated at the stencil 
center and three intermediate points, with one outgoing Riemann variable calculated at a 
stencil edge, and with boundary data for u and p, then this boimdary treatment is imstable. 
A stable boundary treatment is obtained if a nonstandard interpolation is introduced on an 
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eight point stencil with four points on the interior side of the stencil center, and three points 
on the exterior side. The extra degree of freedom of information that is required for eighth 
order accuracy is obtained by adding spatial derivative variables for both u and p on the 
boundaries. The method of undetermined coefficients easily accomodates this information by 
using the derivatives of (6) evaluated at the boimdary in addition to the values of (6) on the 
eight stencil points. For example, in an eight point stencil at the left edge of the conputational 
domain, the solution values on the stencil must be suplemented by the derivative data Ux "_3 
Px ?-3 stencil edge. In this case, the method of undetermined coefficients produces 

the estimate 


“1 = - 13230ur_i +4410 u”_ 2 - 1229u;‘_3 

+1225ur + 11025u;Vi - 2646ur+2 + 490u,\3 - 45ur+4), 

with similar estimates for the other Ua and pa- On this stencil, both the outgoing Riemann 
variable and its derivative are calculated at the boundary point, and two degrees of freedom 
of boundary information must be supplied. The derivative of the outgoing Riemann variable 
is obtained by combining the spatial derivatives of (10). An eight point stencil at the right 
edge of the computational domain is handled in a similar manner. If data is provided for one 
of the variables and its spatial derivative, then this eight point stencil boundary treatment is 
stable and eighth order accurate in both space and time. 

Numerical results are presented in Table C for the intial boundary value problem with 
initial data u(x,0) = 0 zind p(a:,0) = Sin(Tri), for —1 < a: < 1, and with boimdary data 

“ 0) “ Sin(7r(— 1 + <))) and p(l,t) = |(Sin(7r(l — t)) — Sin(7r(l + <))), 
for 0 < <. This problem is the periodic initial value problem from Section III-B in a truncated 
domain with the solutions for u(— l,t) and p(l,t) given at opposite boundaries of the domain. 
This type of boimdaxy data specification is common in computational fluid dynamics. The 
boundary algorithms axe as specified in this section. 


Table C: Data From The Four Algorithms For The ID LEE 
BC - Outgoing Riemann Variable, u(— l,t), and p(l,t) 
u(x,0) = 0 and p(x,0) = Sin(Trx), A = 0.8, M = 0 
Maximum Error in u or p at < = 10 


2 

h 

«10 

Or = 2 

8 

50 

1.870D-01 

16 

100 

4.569D-02 

32 

200 

1.323D-02 

64 

400 

3.500D-03 

128 

800 

8.937D-04 

256 

1600 

2.254D-04 

512 

3200 

5.657D-05 

1024 

6400 

1.417D-05 


Or = 4 Or = 6 Or = 8 

1.170D-02 1.012D-02 8.653D-05 
9.980D-04 5.260D-05 8.499D-07 
8.109D-05 8.786D-07 4.690D-09 
5.540D-06 1.290D-08 2.087D-11 
3.580D-07 1.894D-10 3.752D-13 
2.270D-08 3.346D-12 4.910D-13 
1.430D-09 1.889D-12 
9.115D-11 
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The data in Table C confirms that the propagation algorithms with the bo^dary alg^ 
rithms have the same order of accuracy as the propagation algorithms with periodic bound- 
aries This can also be confirmed by truncation error analysis. The parameters used for the 
computations in Tables A and C are the same, so that the results can be directly compared^ 
The overall error data in Table C is slightly smaller than in Table A, by method and by ^d 
size, most notably for the eighth order method, particularly at = 8. At ^ - 8, the eigh 
order method has a nine point domain, and uses only three stencils for the entue calculation. 
The entire nine point domain is used for computing the values at the central gnd point the 
leftmost eight points are used to calculate everything to the left of the central point, and the 
rightmost eight points everything to the right. Specification of Dirichlet data ^d computation 
of the outgoing Riemann variables provides slightly more accurate results throughout a e 
C than comparable results in Table A with periodic boundary conditions. The computations 
for this initial boundary value problem can be carried out to longer times. For example, in 
computing out to t = 1000 with the eighth order method, if f = 8 is used, then the maximum 
absolute error is 5.405D-05, while if f = 64 is used, then the maximum absolute eiror is 
6 522D-12 These boundary algorithms appear stable, and they clearly have the same order ot 
accuracy as the propagation algorithms, from second to eighth. The bounda^ algorithms c^ 
be viewed as using interior differencing at and near the boundaries. From this viewpoint the 
stencils become heavily weighted towards the interior, much more so than appears possible 
with conventional upwind algorithms. A central feature of these boundary algorithms is the 
simultaneous and consistent calculation of the correct wave dynamics of the solution over an 
interval next to the boimdary. 
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IV: Numerical Algorithms For The Convection Equation In Two Dimensions 


The simplest hyperbolic problem in two space dimensions is the scalar first order linear 
wave equation for u{x,y,t), 

du ^ ^ du du ^ ^ 


dt 


where and My are the constant mean convection velocities in the x and y directions, 
respectively. The initial value problem u(x,y,0) = ui{x,y) for (x, y) e §? x 3ff has the general 
solution u(x,y,t) = ui(x — Mxt,y — Myi) for (x,y) E 3? x 3? and 0 < t. This solution form 
is developed by the method of characteristics, it incorporates the geometry of the solution’s 
behavior, and it can be used to develop exact propagator algorithms for equation (12). In one 
space dimension with a common stencil the exact propagator and Taylor series approaches to 
algorithm development lead to the same algorithm, but for two dimensional problems with 
a common stencil they lead to algorithms that axe distinctly different. Both approaches are 
used and compared in this section with examples of second, fourth and sixth order algorithms. 
A significant new issue for two dimensional problems is the choice of stencil and interpolant. 


IV-A: A General Algorithm Development 


A uniform mesh is used with grid points (xj,yj) = (ih,jh) for integers i and j, where 
h = Ax = Ay is the uniform mesh spacing in space, and with discrete times = nk for 
integer n, where k = At is the uniform time step size. The numerical solution at the mesh 
point (x,-,yj,t„) is denoted by fa u(xi,yj,tn). In two spane dimensions, a polynomial 
spatial interpolation to u around (xj, yj) at t„ cein be written in local coordinates as 

u(x, + x,yj-|-y,t„)««a(x,y)= ^ u^,px°y^, (13) 

{a,;9}eAS 

where AS is the index set for the expansion, and where the coefficients are not indexed with 
respect to the mesh point. All of the undetermined coefficients are simultaneously obtained 
by interpolating a data surface to the solution values on a given stencil. The coefficients can 
be interpreted as spatial derivatives, with 

1 d°‘'^^ua 1 ^ 

The general form of the solution to equation (12) and the interpretation of the expansion 
coefficients can be used to develop a truncated Cauchy-Kowaleskaya series expansion that 
locally approximates the exact solution in space and time, with 

u(xi+x,y, +y,t„ + t) wua(x-Mit,y-My<)= ^ «a,/j(a: - M*<)“(y - Mj,<)^. (15) 

{of,yS}eAS 
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This local polynomial approximation of the exact solution to the global problem is an exact 
solution to the local approximate problem which uses the local spatial interpolant as initial 
data, and it incorporates the correct wave propagation dynamics for (12). 

The exact propagator algorithm form is obtained from the method of characteristics by 
using the locally exact solution form (15), with 


=«a(-Mxik,-Mj,fc) = fnu{xi,yj,tr, + k). ( 16 ) 

{a,y9}e>is 

This algorithm form can also be viewed as a time expansion, with powers of k up to the 
maximum possible a A truncated Taylor series expansion in time produces the algorithm 

form 

<r = 

-y=0 ^=0 

Or 7 


Lyr j ^ 


d~'ua 


7=0 

Or 


(7 - m’- dx^-^dy^ 

d'^u 


fi^ij Vj 9 ^n) 


^=0 

7 


(17) 


5 £ (7 - m'- dx-<-^dy^ w- 


= 'y ] — I" ~ yji^n + ^)? 


^ 7 

y=0 ' 

where Or is the order of the expansion. Both algorithm forms can be rewritten as conventional 
single step explicit finite difference methods, with 

{r,i}G/S 


where IS is the appropriate index set for the stencil that is being used, and where the constant 
coefficients Cr,s 3^® polynomials in the convection velocities Mx and My, and in the grid ratio 
_X = For any given order of accuracy Or, if the expansion coefficients of (13) from a 
fixed stencil are used in both the characteristic based (16) and the Taylor series time expansion 
(17) algorithms, then both algorithms will have the same explicit finite diference index set 
IS, and they both require the same number of floating point operations per grid point per 
time step. In general, the time expansion form of the characteristic based algorithm will have 
time terms that are higher than Or, and the finite different coefficients Cr,a for the two types 
of algorithm will differ. 

The order of the algorithms of either type will depend upon the stencil and interpolant 
that are used for local approximation of the known data surface. Figure 1 presents information 
for stencils and interpolants that can be used for second, fourth and sixth order algorithms in 
two space dimensions. In Figure 1, information for second, fourth and sixth order algorithms 
is presented in the lefthand, center, and righthand columns, respectively. Stencil schematics 



are presented in the top row. The stencils axe all symmetric in space, with a stencil choice 
for fourth and sixth order algorithms. The points at the stencil centers are marked with a 

other necessary points are marked with an ‘o’, and optional points are marked with an 
‘x’. The only second order stencil is the familiar 3x3 square. The two fourth order stencils 
are the 5x5 square, with twenty five points, and this square minus its four corner points, 
with twenty one points. The three sixth order stencils axe the 7x7 square, with forty nine 
points, this square minus its four corner points, with forty five points, and this square minus 
the three points in each comer that are either in the comer or next to it on an edge, with 
thirty seven points. Spatial interpolation coefficient index sets are presented in the second and 
third rows of Figure 1, and represent the index set AS that is used in (13). Each index pair 
represents an expansion coefficient that can be interpreted as a mixed partial derivative in x 
and y, as in (14). The expansion coefficients in the second row are for use with the complete 
square stencil in each case, and the expansion coefficients in the third row are for use with 
the smallest stencil in each case. The index sets in the second and third rows of Figure 1 are 
used in (16) for exact propagator algorithms on either the full stencil or the smallest stencil 
for each order. The fourth row of Figure 1 presents the index sets that are retained in the 
Taylor series time expansion algorithms (17). The coefficients in the fourth row for the Taylor 
series algorithms can be obtained on any of the symmetric stencils that are possible for each 
order, or by other methods. The algorithms in this section axe completely specified by the 
stencils and expansion coefficient index sets in Figure 1, and the general algorithm forms (16) 
and (17). 

It can be shown that if the data for a simulation is independent from one of the coordinate 
variables, and if the mean convection velocity is perpendicular to that coordinate axis, then 
for a given order, algorithm forms (16) and (17) both reduce to the same one dimensional 
algorithm along the mean convection direction. Under these assumptions, the second order 
methods both reduce to the Lax-Wendroff scheme in one dimension, and the fourth order 
methods reduce to the algorithm in Section II. In multiple dimensions with general data and 
symmetric interpolation, algorithms for equation (12) that are derived as tnmcated Taylor 
series time expansions cannot be interpreted as projection backwards along a characteristic 
to its point of intersection with the local data surface at time tn- In effect, Taylor series time 
expansion algorithms introduce an error in the time evolution of the local interpolated data 
surface. 


IV-B: Second Order Algorithms 


The second order algorithms both use the 3x3 square stencil on the left of the top row 
in Figure 1. The expansion coefficient index set for the exact propagator algorithm is given 
on the left of the second row in Figure 1, and the expansion coefficient forms are given in 
Appendix A. The local biquadratic spatial interpolant for u near (x,*, y^) at tn can be written 
as 

2 

u(xi + x,yj + y,tn)!=iua(x,y)= Ua,0X°y^. (18) 

Q,fi=0 
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With this interpolant, the locally exaet algorithm (16) is 

=uoo - k{uioMx + uoiMy) + k'^{u2oMx^ 4- unM^My + uo2Mj,^) 

— k^{U2\Mx^ My + Ul2MxMy“^) + k*(U22Mx My ). 


( 19 ) 


Notice that algorithm (19) has time expansion terms up to k*. The truncation error for 
algorithm (19) is 


1 11 5u , 5tx ,, 


+ 0[h% 


Algorithm (19) is consistent with equation (12), it is second order accurate in space and time, 
3 Xid it is dispersive with no cross derivative terms in its truncation error. 

The truncated Taylor series time expansion algorithm (17) for this stencil is 

= Uoo “ k(uioMx + "I" k^(u2oMx^ + UiiMx^y + '^02^y )j (20) 

which is just the second order Lax-WendrofF method for (12). Notice that algorithm (20) 
includes the terms of algorithm (19) up through P, but does not have P and P terms. 
Algorithm (20) requires the cross difference term «u, and the simplest approximation of uu 
with a symmetric stencil requires the use of the four comer points in the 3 x 3 nine point 
central stencil. The tnmcation error for algorithm (20) is 


1 du du ^ , du 

-(«(x,-,y,-,t„+i) - Uij ^ dy 






6 

d^U ,^ MxMy^ ^^^ d^U 


dx^dy 


-h\ 


dxdy 


+ 0[p]. 


Algorithm (20) is consistent with equation (12), it is second order accurate in space and time, 
it is dispersive, and it has every possible 0[P] cross derivative term in its truncation error. 


IV-C: Fourth Order Algorithms 

The square 5x5 central stencil can be used for fourth order algorithms, with the biquartic 
interpolation 

4 

u(xj + x,yj + y,tn) ««a(a;,y) = ^ Ua,fiX^y^ = ^ Ua,fiX°y^ . (21) 

{a,/9}€A25 
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The locally exact algorithm form (16) with this interpolant is 


=uoo - k{uioMx + uoiMy) + P(«2oM*^ + ujiM^My + uq2Mj,^) 

- P(U3qMx^ + U2lMx‘^My + Ul2MjMy^ + UosMy^) 

+ Ar^(u4oMi^ + usiMx^My + U22Mj^My^ + uisM^My^ + U(nMy^) 

- k\u^iMx*My + U32Mx^My^ + U2iMx^My^ + Ul^MxMy^) (22) 

+ fc®(u42M/Mj,2 + mzMx^My^ + U2AMx^My*) 

- k\uizMx^My^ + m^Mx^My^) 

+ A:*(u44M/M/). 

Notice that this algorithm has time expansion terms up to fc®. 

The second fourth order algorithm uses the reduced stencil with twenty one grid points 
and the modified biquartic interpolation 


u{xi + x,yj + y,tn)^ua{x,y)= ^ Ua,^x°y^, 

{of,^}6A2i 


(23) 


where A 21 is set of indexes in the center of the third row in Figure 1. The coefiicients for A 21 
jure given in Appendix B. With interpolant (23), the locally exact algorithm form (16) is 

=«oo - k{uxQMx + uoiMy) + k^(u 2 oMx‘^ + unMxMy + uo 2 -^»^) 

- P(u3qMx^ + U2iMx^My + Ui2MxMy^ + «03-Mj,^) 

+ k^{U\QMx^ + U3\Mx^ My + V,22Mx^My^ + U\3MxMy^ + UQiMy^^ (24) 

- fc®(u4iM/Mj, + U32Mx^My^ + U2iMx^My^ + ux^MxMy^) 

+ fc®(u42M/M/ + U24Mx^My*). 


Note that algorithm (24) has time terms up to A:®. 

The leading order truncation error for both algorithms (22) and (24) is 


du 


du 


du 


^(u(x,v„W)-u“t.) = ^ + M.gj + M,^ 


+ 0[h^]. 


Algorithms (22) and (24) have the same 0[h^] truncation errors, except that algorithm (24) 
has an additional mixed {3, 3} cross derivative term. The {3, 3} index pair is not in A 21 . 
Algorithms (22) and (24) jure consistent with equation (12), they are fourth order accurate 
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in space and time, and they have leading order truncation error that are dispersive with no 
cross derivative terms. In the form of conventional finite dilference methods, algorithm (22) 
requires twenty five multiplications and twenty four additions at each grid point for each time 
step, and algorithm (24) requires twenty one multiplications and twenty additions. 

The third fourth order algorithm uses the tnmcated Taylor series in time (17), with 

=uoo — k{uioMx + uoiMy) + k^{u2oMx^ + u\\MxMy + uo2My^) 

— fc^(u3oMi^ + U2iMx^ My + u\2MxMy^ + uozMy^) (25) 

-|- k^{u\QMx^ U^\Mx^ My 4" U22^^x "i” d" '^OiMy ). 


Notice that algorithm (25) has the terms up to k^ which are in both algorithm (22) and (24), 
but that it does not have the higher order terms that are used for the ew:curate evolution of 
their local data interpolants. The coefficients in (25) can be obtained from interpolations (21) 
or (23), or elsewhere. If the twenty one point stencil from algorithm (24) is used, then the 
leading order truncation error for algorithm (25) is 


1 j_i du du , , du 




— 9 ' 1\4 ,^^4 d^u 


,4 d^u 

24 ^ dx^dy 


12 dx^dy'^ 


-h\ 


M.M, 
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dxdy^ 


-h\ 


12 dx'^dy^ 


+ 0[h% 


If the twenty five point stencil from algorithm (22) is used to estimate the coefficients for 
algorithm (25), then the 0[h^] tnmcation error terms remain imchanged, and there is a change 
in only the {3,3} mixed 0[h^] term. Algorithm (25) is consistent with equation (12), it is 
fourth order accurate in space and time, it is dispersive, and it has every possible 0[h*] cross 
derivative term in its truncation error. 


IV-D: Sixth Order Algorithms 

Within the 7x7 square stencil there axe three symmetric options for developing algorithms 
that are sixth order, but only two will be considered. Sixth order algorithms will use either a 
7x7 central stencil and the bisextic interpolation 

6 

u{xi + x,yj + y,tn)~ua{x,y)= ^ Ua,fiX°‘y^ = ^ (26) 

[a,fi}€A49 Q,fi=0 
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or the reduced stencil with thirty seven grid points and the modified sixth order interpolation 
u(xi + X, yj + y, tn) W ua(x, y) = ^ u^px^y^. (27) 

{at,^}€A3T 

The index sets A 49 and >137 are given in the righthand colxunn of the second and third rows 
in Figure 1 , respectively. The locally exact solution form (16) produces the algorithms 


“".r = E , (28) 

{q,^}£Aa 9 

with the forty nine point square stencil and interpolant (26), and 

= E (29) 

(a,fi}eA37 

with the reduced stencil and interpolant (27). The leading order truncation error terms for 
both algorithms (28) and (29) are 


1/ / N n 4 -^^ ^ du du 

j(u(xi,v>,(„+i)-<-J )= + + 


+ - M.n“)(4 - M.^A^)(12 - 

+ - Af,n^)(4 - M,^A^(12 - •M,"A^)0 

+ 0[A’]. 


These algorithms are both consistent with equation (12), and sixth order accurate in space 
and time, and they both have dispersive leading order truncation error with no cross derivative 
terms. Notice that algorithm (28) has terms up to while algorithm (29) has terms up to 
A:®. There appears to be little eicctiracy advantage for algorithm (28) with the larger stencil, 
which has approximately a third more grid points. The thirty seven point stencil does not 
have comers, so that special care has to be exercised when treating boundaries and comers 
for problems with nonperiodic boimdaries. 

The sixth order tmncated Taylor series algorithm (16) has the form 




6 


= E 


k-^d^u. , 


(30) 


with terms up to fc®. If the thirty seven point stencil is used to obtain the coefficients for 
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algorithm (30), then the leading order tnmcation error terms are 




dt 


dx 


dy 


M: 




d^u 


5040 


dx‘ 




720 ’ dx^dy 


d'^u 
720 


-h\ 

-h\ 


\ 'Ton 

240 ^ ' dx^d^ " ' 240 "■■ dx^dy 

d^u ^q,Mi My d u 


^M/MyK , 8 '« , an 

O So 2 “ " V o/in a_2a..5 


dxdy^ 

-‘"-‘y \\6 


144 


dx^dy^ 


)X^ 

144 ’ dx^dy^^ 


+ 0[/i^]. 


Algorithm (30) is consistent with equation (12), it is sixth order accurate in space and time, 
and it has a leading order truncation error term that is dispersive with every possible cross 
derivative term. 


IV-E: Numerical Comparisons 

As a numerical test of the algorithms in this section, consider equation (12) with the 
initial data 

it(i,y,0) = Sin(7rx)Sin(7ry), 

for (x,y) G [—1,1] x [-I?!]? 2nd with periodic boxmdaries. Computations will be done with 
Mj. = 1 = My and A = | = ^, for a sequence of grid sizes, and out to the fixed time t = 10. 
N^ nm prir.a.1 data for this problem is presented in Table D. In Table D, is the number of grid 
points in [-1, 1], or per wavelength, nio is the number of time steps required to compute to 
t = 10 at a given grid resolution with A = ^, and the data is the maximum absolute error 
in tt over the entire grid at t — 10. In the column headings of Table D, Or is the order of 
the method, EXg and TSg represent the second order exact propagator and Taylor series 
algorithms (19) and (20) on the nine point 3x3 stencil, EX 21 and TS 21 represent algorithms 
(24) and (25) on the twenty one point stencil, EX 49 and TS 49 represent algorithms (28) and 
(30) on the forty nine point 7x7 stencil, and EXsj and T S37 represent algorithms (29) and 
(30) on the thirty seven point stencil. Note that the square twenty five point 5x5 stencil is 
not used, and that the sixth order Taylor series algorithm (30) is used with both the thirty 
seven and forty nine point stencils. The second order Taylor series method TS 9 is unstable 
at all grid resolutions with A = 0.8 for Mx — 1 and My = 1, so that A = 0.4 is used to obtain 
it’s data. The large error for the second order Taylor series algorithm at | = 8 is due to 
dispersion errors from poorly resolved data, just as for similar cases among the data in Table 
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A. Calculation of error reduction exponents similar to those conducted for the data in Table 
A confirms the order of accuracy for each of the methods used to produce Table D. 

Table D: Data For The Scalar First Order Wave Ex^uation In 2D 
u(x, y, 0) = Sin(7rx)Sin(7ry), A = 0.8, Mx = 1 = My 
Maximum Error In u at < = 10 




Or = 2 

Or = 2 

Or = 4 

Or = 4 

Or = 6 

Or = 6 

Or = 6 

Or = 6 

2 

h 

nio 

EX^ 

TS^* 

EX 21 

TS21 

EX49 

TS49 

EXz7 

TSz7 

4 

25 

0.999 

0.959 

0.978 

0.675 

0.670 


0.733 

0.841 

8 

50 

0.895 

1.414 

0.156 

0.365 

1.8D-02 

2.0D-02 

2.2D-02 

2.2D-02 

16 

100 

0.322 

0.387 

9.6D-03 

4.2D-02 

2.8D-04 

2.3D-04 

3.1D-04 

2.6D-04 

32 

200 

8.0D-02 

8.6D-02 

5.4D-04 

2.8D-03 

4.1D-06 

2.9D-06 

4.3D-06 

3.2D-06 

64 

400 

1.9D-02 

2.0D-02 

3.2D-05 

1.8D-04 

6.0D-08 

4.0D-08 

6.2D-08 

4.2D-08 

128 

800 

4.7D-03 

4.7D-03 

1.9D-06 

l.lD-05 

9.1D-10 

6.0D-10 

9.2D-10 

6.1D-10 

256 

1600 

1.2D-03 

1.2D-03 

1.2D-07 

7.0D-07 

1.4D-11 

2.8D-06 

1.4D-11 

9.2D-12 

512 

3200 

2.9D-04 

2.9D-04 

7.3D-09 

4.4D-08 

2.2D-13 

2.9D-05 

2.2D-13 

4.0D-11 


* Note that A = 0.4 for the data in this column. 


The sixth order Taylor series methods become inaccurate at fine grid resolutions with 
A = 0.8, neeir ^ = 256 for TS 49 , near j = 512 for TSzr. If the grid ratio is reduced, then 
these Taylor series methods produce accurate results at these grid resolutions. The fine grid 
errors with the T ^49 and T S37 methods appear to be due to the excitation of parasitic steady 
state periodic solutions of equation (12). Notice that the Taylor series algorithm on the larger 
square stencil shows signs of inaccuracy at a coarser grid resolution than the algorithm on 
the reduced stencil. Numerical experiments show that exact propagator methods are stable if 
both one dimensional CFL numbers are less than one in absolute value, but that the Taylor 
series methods have more restrictive stability constraints. At eaxdi level of grid refinement, 
the data in Table D shows errors from the exact propagator sind Taylor series methods on the 
same stencil that differ by at most a factor of about 5. The largest differences are in the fourth 
order results, with smaller errors from the exact propagator method. Recall that both types 
of method can be recast as conventional finite difference algorithms, and that in this form the 
munber of operations depends only upon the stencil. The exact propagator methods appear 
to be at least comparable in accuracy to the Taylor series methods, if not more accurate, and 
they appear to be more robustly stable, without requiring more operations. The choice of 
largest or smzillest possible symmetric stencil does not appear to be significant for accuracy. 

The data in Table A for the one dimensional linearized Euler equations, and the data in 
Table D for the two dimensional convection eqtxation both come from periodic initial value 
problems with pure frequency sine wave initial data at the same wavelength. The two dimen- 
sion2il convection calculations are run with the velocity (Mx,My) = (1, 1), which is at a 45® 
angle with the grid lines. The close agreement between the error levels of the comparable one 
and two dimensional results suggests that convection skew to the grid does not introduce sig- 
nificant errors if these algorithms are used. Note that both algorithm types use time evolution 
that reflects the genuinely multidimensional nature of the convection equation. The Taylor 
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series in time methods include the cross derivative terms required for higher than first order 
accuracy in time, and the exact propagator methods directly incorporate the characteristic 
behaviour of the solution. Both method types use genuinely multidimensional interpolation 
for the local approximation of the solution surface in order to obtain the expansion coefficients 


that are needed for accurate time evolution. 

K interpolation is used with the data in Table D, then in order to achieve a maximim 
absolute error of less than 5.0 x 10"^ at t = 10, the second order methods require ^ > 216, the 
exact propagator fourth order method requires ^ > 17, the Taylor senes fourth order method 
requires ^ > 30, and the sixth order methods require ^ > 8. For a grid ratio of A = ;^, the 
number of time steps nio required to compute to t = 10 is nio = Consequently, the tot^ 
munber of multiplications required for each computation on the doinain (x, y) € [-1, Ij x [-1, IJ 
with j TTi p«b points per unit interval, using a stencil with n, grid points, and for nio time 

steps is 


total multiplications = nio^h^ris = 


40h^n. 


In order to meet the 5.0 x 10"^ error bound at t = 10, the ratio of the number of multiplications 
required by the different order methods is approximately 


40 X 216^ X 9 
40 X 173 X 21 


879.1, 


for either of the second order methods to the exact fourth order method, and 


40 X 17^ X 21 ^ . 

40x83x37 - “ 

for the exact fourth order method to either of the sixth order methods. If the Taylor senes 
fourth order method is used, the ratio of the number of multiplications required by the different 
order methods becomes approximately 160.0 for the second order methods to the fourth order 
method, and approximately 29.9 for the fourth order method to the sixth order methods. 
These comparisons show that the efficiency advantage of higher order methods is increased for 
increased spatial dimension, as well as for increased absolute simulation time and decreased 
maximum error Umit. The greatest payoff by fax is in the relative efficiency shown by the 
exact propagator fourth order method when compared to the second order methods. 
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V: Numerical Algorithms For Linearized Euler Ekiuations In Two Dimensions 


The linearized Euler equations in two space dimensions can be written as 


du du du dp 

^ ^ = 0 , 

dt ‘dx ^dy dx 

dv dv dv dv 


dt 


dx 


dp dp dp du dv 


(31) 


where p is the pressure, (u,u) the distmbance velocity, and (Mx,My) the constant mean 
convection velocity. This formulation is for the isentropic case, and is nondimensionalised in 
terms of the Mach number. The linearized Euler equations aire essentially multidimensional 
because they cannot be diagonalized and transformed into a simpler set of decoupled equations, 
and wave propagation is along characteristic surfaces instead of characteristic cmves. 

A imiform mesh is used, with mesh sizes h and k in space and time, and with numerical 
solutions denoted by u”j, v^j, and p"y. Polynomial spatial interpolations to u, u, and p are 
written in local coordinates around (xi,yj) at with the form (13). Second, fourth and 
sixth order methods will use the nine, twenty one, and thirty seven point stencils in Figure 1 . 
Interpolation coefficients are obtained by the method of tmdetermined coefficients, and can 
be interpreted in terms of spatial derivatives as in (14). Exact polynomial solution forms for 
the linearized Euler equations can be derived by substituting the expeinsion forms 


Or 

u{xi+x,yj + y,tn + t)^ua{x,y,t)= ^ ^Ua,fi,yX°‘y^t'^ , 

{a,P}£AS-r=0 

Or 

v(xi + X,yj + y,tn + t) ^ va(x,y,t) = ^ Y^Va,0,yX°y^f^ , (32) 

{a,P}£AS 1=0 
Or 

p{xi+x,yj + y,tn + t)fnpa(x,y,t)= , 

{a,p}£AS 1=0 


into system (31), and obtaining all the terms with 7 ^ 0 by reqmring system (31) to be 
satisfied for all 2 ;, y and t. Coefficients with 7^0 are equivalent to time derivatives, and 
the resulting polynomial solutions are expressed entirely in terms of the spatial expansion 
coefficients. The exact polynomial solution forms with the local spatial interpolants as initial 
data give exzict propagator algorithms, and automatically incorporate into the solutions the 
correct local multidimensional wave propagation dynamics for the local spatial interpolants. 


V-A: Second Order Algorithms 


The second order methods approximate u, v and p with the biquadratic spatial interpola- 
tion (18) on the 3x3 stencil. Formulas for the u expansion coefficients are given in Appendix 
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A. The second order exact propagator method has Or = 4 in (32), with the solutions given in 
Appendix C. The exact propagator algorithm is 


=ua(0,0,fc) 


=«oo 

+ k{—pio — MyUQi — M*txio) 

+ k^{MyPii + 2M^P2o + My^U02 + MxMyUxi + (1 + M^^)u20 + 

+ + -^»^)Pl2 - 2 MxMj,P21 

- M:cMy^Ui2 - {My + Mx^Mj,)u21 - MyVi2 - MrV2l) 

+ + 2M^My^)p22 + ( J + M^^My^)u22 + 2M^MyV22), 


(33a) 


=va{0,0,k) 

=Voo 

+ k { — poi — MyVQi MxV\q') 

+ k^{2MyP02 + MxPii + ^Uii + (1 + MJ)uo2 + MxMyVii + M^V2o) 

+ k^{—2MxMyPi2 ■" (3 + MI)p2i 

— MyUi2 — MxU2\ — (Ml + MxMy)vi2 — MlMyV2i) 

+ k*{{\My + 2MlMy)p22 + 2MxMyU22 + (^ + M^ + M^MI)v22), 


and 


Pi,V 


(336) 


=Poo 

+ k{—MyPoi — MxPio ~ ^10 ~ t^oi) 

+ fc^((l + My)po2 + MxMyPii + (1 + M|)p2o 

+ MyUii + 2 MxU 2 o + 2 MyVQ 2 + MxVii) (33c) 

+ P{-{Mx + MxM2)pi2 - (My + MlMy)p2i 

-i\+ MJ)ui2 - 2MxMyU2i - 2MxMyVi2 - (^ + MI)v2i) 

+ k*{{^ + M^ + Mj + MlMl)p22 

+ (2Mx/3 + 2MxMI)u22 + (2Mj,/3 + 2MlMy)v22). 

Notice that algorithm (33) has the form of a time expansion with terms up to If algorithm 
(33) is written in the form of a conventional explicit finite difference method, then it requires 
twenty seven constant coefficients for each of the three equations. The time expansion form 
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is easier to develop and program, and is more flexible, but the finite difference form is more 
computationally efficient. 

The tnmcation errors for algorithm (33) are 


_ (3 + 


dx 
+ 0[h% 


(34a) 




+ '>"0(g)(l-(l + 3W,^)A^) 

+ 0[h% 


i(p(xi,yj, t„+i) - ?"+■) = + 

+ A^0(^)(1-(3 + MJ)A^) 

+ fc^0(i)(l-(l + 3M|)A=) 

+ '>^0(i)(l-(l + 3M,^)A») 
+ 0[4’|. 


(346) 


(34c) 


Note that the tnmcation errors (34a) and (34c) for u and p are the lowest order errors of 
the second order method from Section III plus one additional third derivative term eaeh. The 
truncation errors (34) clearly show that algorithm (33) is second order accurate and dispersive. 

The second order Taylor series time expansion algorithm on the 3x3 central stencil can 
be obtained by taking the time terms up through k“^ in (33). If symmetric crossdifferencing is 
used to obtain un, Un and pn, then the second order Taylor series time expansion algorithm 
requires the same number of floating point operations per grid point per time step as algorithm 
(33). With this differencing, the truncation error for u from the second order Taylor series 
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time expansion algorithm is 


w.) - <r) 


,d^p,M, 


+ +MJ)M,A=' + h-‘-^(=^)M,Ml\-‘ 


dx^dy 2 


+ h 


dx^dy'^ 2 
d^p 


dxdy 




+ 0[h% 


This tnmcation error for u from the time expansion method contains the error (34a) for u 
from the exact propagator zilgorithm (33) plus an additional six cross derivative terms. The 
tnmcation errors for v and p have a similar complexity with cross derivative terms that are 
absent from (34b) or (34c). Note that algorithm (33) correctly incorporates the local wave 
dynamics for its spatial interpolant, but that the Taylor expansion in time does not. 


V-B: Higher Order Algorithms 

For the sake of brevity, only an outline will be given of higher order methods, but with 
sufficient detail to ensme their specification. The local approximation of u, v ^d p is done 
for the fourth order algorithms with the twenty one point stencil and interpolation (23), and 
for the sixth order algorithms with the thirty seven point stencil and interpolation (27). With 
these approximations the fourth order exact propagator method has Or = 6 in (32), and the 
sixth order has Or = 8. The exact solution forms are used to obtain the numerical algorithms 
with = ua(0,0,ib), = va(0,0,k), and ^ pa(0,0,fc). The time expansion 

forms of the fourth order algorithm are given in Appendix D. Notice in Appendix D that the 
algorithms for and are symmetric in u and v, if the role of z and y are interchanged. 
The tnmcation errors for the fourth order exact propagator algorithm axe 


^(u(xi,yj,i„+i) 


u"+i) = - 5X^(3 + M") + A^(5 + lOM^ + M^)) 

120 ox^ 

__^^Mj,(4 - 5X^M^ + A^MJ) 

120 ay® ^ 

- 5X^(3 + M^) + A"(l + 10M2 + 5M^)) 

120 dx 


+0[/i®]. 


(35a) 
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(356) 





1 

k 


J/jf? ^n+l) PiJ ) 


-^0Af,(4 - 5A^(3 + Ml) + A‘(6 + lOM," + Mj)) 

-^5(4 - 5A^(3 + Ml) + A^l + lOMl + 5Mj)) 

+0[A'l, 

- 5A^{3 + Ml) + A<(6 + lOMl + M^) 

~Tm ” 5A^(3 + MI)+ A^(1 + lOMl + 5Ml)) 
-^0M,(4 - 5A^(3 + *(^) + A<(5 + lOMj + Mj)) <36c) 

-^1^(4 - 5A“(3 + Ml) + A*(l + lOM,^ + SJW*)) 

+0|A=]. 


Note that the truncation errors (35a) and (35c) for u and p in two dimensions are the lowest 
order error terms of the fourth order method from Section III plus an additional fifth derivative 
term. The truncation errors for both methods are dispersive and completely lacking in cross 
derivatives, and they show that the methods are fourth and sixth order accurate in space and 
time, respectively. 

The fourth order Taylor series time expansion algorithm can be obtained by dropping 
the and terms from the exact propagator algorithm in Appendix D. Similarly, the 
sixth order Taylor series method can be obtained by dropping the k^ and A:® terms from the 
sixth order exact propagator algorithm. The truncation errors for these methods confirm 
their order of accuracy in both space and time. None of the Taylor series algorithms can be 
interpreted as correctly incorporating the wave propagation dynamics of their interpolants, 
and they exhibit a plethora of cross derivative terms in their lowest order truncation errors 
because these derivatives from the spatial interpolation axe not incorporated in the local time 
evolution. The number of error terms that the Taylor series algorithms add to the lowest 
order truncation errors from the exact propagator algorithms increases with the order of the 
aJgorithm. 


V-C: Numerical Comparisons 


As a numerical test of the algorithms that axe discussed in this section, consider system 
(31) with the initial data 

p(x,y,0) = Sin(7rx)Sin(7ry), 
ii(x,y,0) = 0, 
u(x,y,0) = 0, 
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for (x,y) G [-1,1] x [-1, 1], with periodic boiindaxies. The exact solution is 

p{x,y,t) = Cos(v^xt)Sin(7r(x - Mxf))Sin(7r(y - 

u{x,y,t) = ^Sin(\/27rt)Cos(7r(x - Mxt))Sin(7r(y - MJ)), 

V2 

v{x,y,t) = 4=Sin(\/27ri)Sin(7r(x - Mxt))Cos{n{y - 

V2 

for (x,y) e [-1,1] X [-1,1] and t > 0. This problem is the direct two dimensional extension 
of the problem used in Section III-B for the numerical comparison of algorithms for the one 
linearized Euler equations, and similar data is used for the numencal tests in 
Section IV-D. The algorithms in Section III for the linearized Euler equations in one space 
dimension have the stability constraint A = f < and in that section A = 0.8 is used with 

mean convection velocity M = 0. Two dimensional analogs of the one dimensional stability 
constraint include replacing 1 + |M| by either 1 + (M^ + or by 1 + Max{|Mx|, lMy|}. 

For Mx = 1 = My, this becomes either A < 1/(1 + •v/2) 0.414, or A < 1/(1 + 1) = 0.5. With 

Mx = 1 = My, the exact propagator algorithms for equation (31) are stable with A = 0.5, the 
second order Taylor series algorithm is not stable with A = 0.4, but it is stable with A = 0.25, 
and the fourth and sixth order Taylor series algorithms are not stable with A = 0.5, but they 
are stable with A = 0.4. The computations in this section will be on a sequence of grid sizes 
out to the fixed time t = 10, generally with A = J = 0.4, except with A = 0.2 for the second 
order Taylor series algorithm. Numerical data for this problem is presented in Table E. 


Table E: Data For The Linearized Euler Equations In 2D 


p{x,y, 0) 

= Sin(' 

7rx)Sin(7ry), 

u(x,y,0) 

= 0 = u(a 

^y,o), A 

= 0.4,Mx 

= 1 = M, 


Maximum Error In p at 

t = 10 





Or = 2 

Or = 2 

Or = 4 

Or = 4 

Or = 6 

Or = 6 

2 

h 

«10 

EX^ 

T5g* 

EX 21 

TS 21 

EX^-r 

TSsj 

4 

50 

0.885 

0.851 

0.931 

0.917 

0.735 

0.763 

8 

100 

0.888 

1.054 

0.168 

0.262 

2.2D-02 

1.7D-02 

16 

200 

0.344 

0.524 

9.9D-03 

2.5D-02 

3.0D-04 

1.9D-04 

32 

400 

0.102 

0.145 

6.9D-04 

1.7D-03 

5.3D-06 

3.4D-06 

64 

800 

2.7D-02 

4.0D-02 

4.5D-05 

l.lD-04 

8.8D-08 

5.9D-08 

128 

1600 

6.9D-03 

l.OD-02 

2.9D-06 

6.8D-06 

1.4D-09 

9.7D-10 

256 

3200 

1.7D-03 

2.6D-03 

1.8D-07 

4.3D-07 

2.3D-11 

5.1D-12 

512 

6400 

4.4D-04 

6.6D-04 

1.2D-08 

2.7D-08 

1.6D-11 

5.0D-12 

'' Note that A = 

0.2 for the data in this column. 





In Table E, ^ is the number of grid points in the interval [—1,1] or per wavelength, nio 
is the number of time steps reqmred to compute to t = 10 at a given grid resolution with 
A = — , and the data is the maximum absolute error in p over the entire grid at t = 10. In the 
colunm headings of Table E, Or is the order of the method, and EXm and TSm represent the 
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exact propagator and Taylor series algorithms on the m point stencil. Calcvdations of error 
reduction exponents similar to those conducted for the data in Table A confirm the order of 
accuracy for each of the methods used to produce Table E. The error levels in Table E are 
simil£ir to those in Table D when compared by method and by grid size. The large error for 
the second order Taylor series algorithm at ^ = 8 is due to dispersion errors from poorly 
resolved data, just as for similar cases among the data in Tables A and C. If these methods 
are compared by order, then the exact propagator and Taylor series methods are similaa* both 
in the computational effort that they require, and in the errors that they produce, but the 
ex 2 ict propagator methods have less severe stability constraints. For these computations the 
total number of multiplications can be calculated as in Section IV, with 


total multiplications = nio4h^ngnen„ 


A 


where this problem involves ng = 3 equations using data from n„ = 3 variables in each 
equation. Interpolation of the data from the exact propagator methods in Table E implies 
that in order to achieve a maximum allowable error of 5.0 x 10“"* at t = 10, a grid resolution 
of ^ = 248 is required for the second order method, ^ = 20 for the fourth order method, and 
^ = 8 for the sixth order method. The ratio of the total number of multiplictions required 
to meet this error constraint at this time is approximately 817.1 for the second to the fourth 
order methods, and 8.9 for the foiirth to the sixth order methods. These results are similar to 
the comparisons obtained in Section IV, showing greater efficiency for higher order methods. 


V-D: Numerical Computations With Boundaries 


Algorithms require boundary conditions in order to be useful. This subsection presents 
the results of computations with both real and artifical outflow boundaries using the fourth 
order exact propagator algorithm. Details on boundary treatments are given in [15] and [12]. 
Consider the hnearized Euler equations (31) with the initial data 

p(r,y,0) = exp[-/n(2)( - )], 

u(z,t/,0) = 0 = «(x,t/,0), 

with Afx = 0.5 aind My = 0, and on the computational domain (x, y) G [—100, 100] x [0, 200], 
where there is a wall at y = 0, a computational inflow boimdary at x = —100, a computational 
outflow boundary at x = +100, and a far field boundary at y = 200. This problem represents 
a Gaussian pressure pulse near a wall in a Mach 0.5 flow parallel to the wall. Calculations 
will be done with h = l and k = 0.25 out to < = 150. 

The wall boimdary conditions are 


u = 0, 
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Xhe normal derivative of p is evaluated on the wall with fourth order interior differencing 
on a one sided five point stencil. In a wall boimdary cell, the fourth order method is u^d 
to produces values for the three variables at the cell center and at one intermediate point 
between the cell center and the wall, and for u on the wall. All of these solution values 
are foimd in essentially the same way as for the one dimensional equations in Section III-C, 
using the algorithm in the form of a Cauchy-Kowaleskaya expanion in space and time, with 
expansion coefiicients that are obtained from one interpolation on the boundary cell. The 
form of the algorithm as an expansion in space and time is evaluated at (0, 0, k) to obtain the 
three solution values for the cell center, at (0, k) for the intermediate point values, and at 
(0, — 2/», k) for the value of u on the wall. 

New outflow conditions that have been derived by Hagstrom [15] are used at the outflow 
boundary x = +100. If system (31) is diagonalized in the x direction, or normal to the 
botmdary, then the primitive variables p, u and t; can be replaced with 


ri=u-p. 


T2 = V, 


V3 = U+p. 


For subsonic flows, r\ is conventionedly viewed as coming into the computational domain from 
the +x direction, and both T 2 and ra as going out. The variable ri is obtained on the outflow 
boundary bs the solution of the system 


— Mx ^ (/i + /z + + 32) — 

at ay 

^ = -^-(1 - 
dgk 


dy 

dgk 0k , 


2 ’ 




dy 

d^p 


dt 


dy 


dy 


2 ’ 


for jfc = 2, where ai = Cos(f), 02 = Cos(^), = |Sin(f), and 02 = |Sin(^). Notice 

that this boimdary condition requires no geometric information about the disturbance, neither 
globally nor locally. Details are given in [15]. The algorithmic implementation of this boundary 
condition is similar to the algorithm for the wall conditions. In an outflow boundary cell, the 
sj>ace and time expansion is evaluated at (0, 0, k) for the three solution values that are needed 
at the cell center, at (h, 0, ifc) for the three values at the intermediate point, and in combination 
at (2h,0,Jb) for the values of T 2 and that are needed on the outflow. Separate expansion 
forms are developed on the boimdary for ri, fi, f 2 , gi, and g 2 , using the expansion forms for p 
and V. These expansion forms use a stencil of width five, and are designed to be fourth order 
accurate in both space and time. The system for ri and the fk and gk is one dimensional on 
the boundary line, forced by the evolution of p and v. The values of ri and the fk and gk 
are all initialized at 0. At the intersection of the wall y = 0 and the outflow x = +100, fk is 
obtained by solving its PDF using interior differencing, while gk is obtained from 


dgk 

dy 


oy 


(36) 
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for A: = 1 and 2. At the far field y = +200 on the outflow x = +100, fk=0 and gt is obtained 
by solving its PDE using interior differencing, for k = 1 zind 2. Details are given in [12]. 

Chatfacteristic boundary conditions are used on the inflow boundary at x = —100, and 
on the fair field boundary at y = +200, with outgoing Riemann variables solved using interior 
differencing and the appropriate momentum equation, zind incoming Riemann variables set 
equal to 0. For these two boundaries, the Riemann variables are determined by a one dimen- 
sional diagonahzation of the system in the direction normal to the boimdary. The propagation 
edgorithm, the wall conditions, and the outflow boimdary algorithm are all implemented as 
fourth order methods, in both space and time, with central stencils having a stencil width of 
five grid points, unless they are offset by interior differencing at a boundary. The inflow and 
far field boundary algorithms are implemented with a second order method using a stencil 
with a three point width, since nothing ever happens at these two boundaries for t < 150. 

The results for this wall pulse problem are presented in Figure 2, which shows pressure 
contours at t = 15, 45, 75, and 150. These results have been briefly described in [11]. In 
Figure 2a at t = 15 the expanding pressure wave just begins to touch the weJl. In Figure 
2b at t = 45 there is a very evident reflecting wave expanding behind the expansion front, 
with interferance patterns where they interact near the wall. The expanding wave reaches 
the artifical outflow boundry at about t = 60, and shows no disturbance at this moment 
of first contact. In Figure 2c at t = 75 the two waves have already passed through the 
artifical outflow, and have progressed from being peirallel to the outflow boimdary at their 
first cont8M;t, to being at approximately 45° with the boundary. Notice that there are no 
evident disturbances in the wave front contours near the boundary at t = 75, even though the 
expanding waves have just passed through the intersection between the wall and the outflow. 
This type of intersection can create a transient error that is propagated with the waves as a 
reflection treiiling backwards from the intersection between the solution and the boundary. In 
Figure 2d at t = 150 the two wave fronts have become nearly perpendicular to the outflow 
boundary, and the pulse center has convected to (x,y) = (75,25). The solution downstream 
from the artificial outflow boundary at x = 100 continuously effects the analytic solution for 
X < 100 over the time interval 60 < t < 150 during which there is significant downstream 
solution structure. The plot sequence shows a correct evolution on the computational domain, 
with solutions that are symmetric in x and have round wavefronts, with no visible perturbation 
of the solution right up to the boundary, and with no evident reflection from the boundary 
or the wall, either as a transient or cumulative effect. This data reflects both the accuracy 
with which the propagation algorithm simulates the fully multidimensional wave dynamics of 
the Hnearized Euler equations, and the unobtrusive quality of the artifical outflow boundary 
condition. 

It could be eirgued that the structure of the wall pulse solution is so simple near and 
downstream from the outflow boundary that the quality of the outflow boundary is not seri- 
ously tested. A more severe test of the algorithm and outflow condition is given by placing 
the initi 2 d Gaussieui pressure pulse in a duct. Consider the same initial data for the Gaussian 
pressure pulse but on the computational domain (x,y) € [—100, 100] x [0,50], with two walls 
at y = 0 and y = 50, with an inflow at x = —100, and an outflow at x = +100. The fourth 
order exact propagator algorithm is used, with h = 1 and k = 0.25 on a uniform grid. The 
only chjinge in the outflow boundary is that at the intersection of the wall y = 50 and the 
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outflow X = +100, 9k is obtained by solving its PDE using interior differencing, while /t is 
obtained from (36) for ifc = 1 and 2. The results for this duct pulse problem axe presented in 
Figure 3, which shows pressure contours at t = 15, 45, 75, and 150. Notice that the wave front 
hits the two walls and is reflected back and forth several times. Each time this happens, a 
more complex structure of interference patterns is visible in the pressure contours. The same 
story is visible in these calculations as from the wall pulse problem, but now the pressure 
structures axe very much more complex, both as they axe successfully passed through the 
numerical boimdaxy, and as they continue to effect the evolution of the pressure structiires 
within the numerical domain because of their virtual representation by means of the bound- 
ary condition. In Figure 3a at t = 15 the expanding pressure wavefront begins to touch both 
walls. In Figure 3b at t = 45 there axe two reflecting wavefronts that have almost reached 
each other in the center of the duct. In Figure 3c at t = 75 the two reflecting wavefronts 
have passed through each other while propagating completely across the duct, and have been 
reflected once again to pass through each other a second time and reach the wall from which 
they were originally reflected. The complex interference pattern at the upstream end of the 
disturbance has a rich structure created by the crossing and recrossing of the reflecting waves. 
Notice at t = 75 that approximately half of the complex structure at the downstream end of 
the disturbance has already passed out of the computational domain. In Figure 3d at t = 150 
the two wavefronts have gone through further reflections and interactions, creating a very 
complex upstream interference pattern. Notice that the entire interference pattern is missing 
from the downstream end of this plot, since it has passed entirely out from the computation^ 
domain. If a transperancy is made of this plot, it can be turned over and the contour hues in 
the center will allign perfectly, with the transparency contours near the center of the domain 
exactly overla 3 ring the plot contours near the boimdaxy. This symmetry in the plot shows that 
even though nearly half of the solution structure is outside of the compuational domain, there 
is no visible effect on the solution inside the domain. The lack of discemable effect from the 
outflow boundary occurs in spite of the fact that the propagating wavefronts in the computa- 
tional dommn are virtually perpendicular to the outflow boimdaxy. This computation shows 
the imobtrusive and robust quality of the artificial outflow boundary condition [15]. 
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VI: Summary And Discussion 


This paper presents two techniques for developing zJgorithms for hyperbolic systems in 
multiple space dimensions, the use of local exact polynomial solutions, or exact propagators, 
and the use of miiltivariate Cauchy- Kowaleskaya expansions, or truncated Taylor series. Both 
techniques of algorithm development can be viewed as approximating the solution of a hy- 
perbohc system as a whole, rather than individual derivative terms in separate equations. 
The multivariate polynomial solution approximations that axe used by these techniques in- 
clude various cross derivative terms, which are needed for high order accuracy, and improved 
isotropy and stability. Both of these two techniques are used to develop algorithms for the hn- 
ear convection equation and the linearized Euler equations, in one and two space dimensions. 
The exphcit single step algorithms developed by these two methods use symmetric stencils 
and axe dispersive, they use data from only one time level, they have the same order of accu- 
racy in both space and time, and they can be extended to arbitrarily high orders of accuracy. 
Algorithm examples axe given with up to eighth order accuracy in one space dimension, and 
sixth in two dimensions. For each order of method, and for each type of algorithm, the choice 
of large or small symmetric stencil in two space dimensions does not appear to be significant 
for accuracy. The higher order methods are more efficient in terms of the total number of 
floating point operations required to compute to a fixed simulation time with a stated error 
bound, and they can require less operations than a low order method by factors of up to sev- 
eral orders of magnitude. The relative efficiency of the higher order methods increases with 
either the simulation time or the dimension of the system domain, and with decreases in the 
error boimd. The exact propagator methods appear to be at least comparable in accuracy to 
the Taylor series methods, if not more acciurate, and they appeeir to be more robustly stable, 
without requiring more operations. Exact propagator algorithms incorporate the correct local 
multidimensional wave propagation dynamics for their polynomial spatial interpolants, and 
they can be viewed as a generalization of the method of chsuacteristics to nondiagonalizable 
hyperbolic systems in multiple space dimensions. In particular, exact local polynomial solu- 
tions are shown for the linearized Euler equations in two space dimensions. Stable high order 
boundauy conditions axe possible using a consistamt calculation of the correct wave dynamics 
of the solution simultaneously from the center to the exterior edge of a grid cell on the bound- 
ary. High order bovmdaxy conditions are developed for the linearized Euler equations in one 
space dimension, and demonstrated in two. 

There are several evident extensions, and issues that are not addressed in this paper. 
Analysis of the algorithms is being done, particularly for stability, with an investigation of the 
different stabihty constraints for the exact propagator and Taylor series methods in two space 
dimensions. The form of the algorithms axe independent from the spatial interpolation that 
is used, so that various methods can be used for estimating the spatial expansion coefficients. 
In particular, an implementation on a nonuniform or unstructured grid is being prepared with 
an appropriate spatial interpolation, and it should not incur any degredation in the order of 
accuraw:y. Implementation in other coordinate systems is being explored. Complete details 
of the outflow boundary algorithm axe being prepared for publication [12]. In particular, 
compatibihty conditions for the juncture of two artifical boundaries are being developed, and 
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implementations for an artifical inflow are being prepared. These two approaches to algorithm 
development are being extended to other flrst order linear systems, such as for Maxwell’s 
equations. The Taylor series method is being extended to variable coefficient cases. 
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Appendix A: Biquadratic Expansion Coefficients For A 3 x 3 Stencil 

The spatial approximation to u around {xi,yj) at t„ is 


u{xi + x,yj + y,t„)!^ua{x,y)= ^ Ua,px°y^, 

Q,p=0 


with difference coefficients 


“00 = 


“11 “ ““H-i.jf-i -“"-1,7+1 +“r-i,i-i)> 

“21 — 4^(“i+i,i+i — «iVi,i-i +2u”j_i — 

“02 = ~ 

“12 = ^^(“”+1,7+1 “ + “"+1,7-1 “ “"-1,7+1 + 2“r-l,7 ~ “^-1,7-1)? 

“22 = ^(“r+1,7 + 1 - 2“"+lJ + «"+l,j-l - 2 “" 7+1 +^“i !7 “ 

o..n t ^.n \ 


+ “"- 1 , 7+1 “ 2u"_ij + 
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Appendix B: Expansion Coefficients For A Modified Biquartic On A Twenty One 
Point Stencil 


The spatial approximation to u around at is 

ua{x, y) =uoo + uiqX + u^qx^ + u^qX^ + 

+(uoi + UuX + U2lX^ + + U4ix‘*)y 

-h(l ^02 + Ui 2 ^ + + ^ 32 ^^ + ^i2^*)y^ 

+(uo 3 + U 13 X + U23^^)y^ 

+(U04 + U143; + U24®^)y^ 

{a,^}GA 2 i 


with difference coefficients 


tXoO — ’ 


= + i6«r+ij - <+2j). 


U20 — 




U 40 = 


2 ^ W-2,j - 4“?-,., + 6»"; - + “?+2,i). 

•*01 = “ 3“i‘.i-i + 3“i!i+i ” “",>+ 2 )’ 

“11 = 24 ^*^ 
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“21 = - “"-2,i+l 
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U23 

«04 

«14 

«24 


^ - {(uu.i-2 - + 2ur-.,i+. - «r-i,>+2 

— 2u"j_ 2 + 4u”j_i — 4u"j^i + 2ui j^2 
+ «"+l,j-2 - 2ur+i,j-l + 2<+ij+i - <+i,j+2)> 

2^(“U-2 “ 

- ((>*?-., i-2 - + 6ur_.,i - 4«”_, j+. + uu„ 

— + 4u”^i j + 4«j^j — «i^.ij+2)» 

^^(*^r-i,j-2 - - 4u”_ij+i + «r- 

— 2u"j_ 2 + — 12u”j + 8uj — 2uj , 4.9 

+ j_2 — 4u".^i j_i + — 4uf.f.i_j,^.] 


•l,J+2 


i,i+2 

,i+2 

r+ij+i + “?+i,i+2)- 
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Appendix C: Exact Second Order Polynomial Solution For Equation (31) 

The locsd exact solution to the linearized Euler equations that uses a biquadratic spatial 
expansion for initial data is given for u by 

ua{x,y,t) = uoo 

+ uio(a: - Mxt) + uoi(y - Myt) 

+ U2o((3^ — Mxi^ + — Myt') + ~ Myt)^ 

+ U 2 i{{x - Mxt)^ + t^){y - Myt) + Ui2(x - Mxt){y - Myt)^ 

+ U 22 {{(x - Mxty + t^){y - MytY + <^/6) 

+ vii{t^ (2) 

+ V2i(x - Mxt)t^ + ui2(y - Myt)t^ 

+ V22{‘^{x-Mxt){y-Myt)e) 

+ Pio(— 0 

+ P2o(-2(x - Mxt)t) + Pii(-(y - Myt)t) 

+ P2i(-2(x - Mxt){y - Myt)t) +Pi2(-(y - Myt)"^t - 1^/3) 

+ P22(-2(x - Mxt){{y - Myt)^ + t^/3)t), 

for V by 

va{x,y,t) = Uoo 

+ uio(x - Mxt) + uoi(y — Myt) 

+ V2o(x - Mxtf + un(x - Mxt)(y - Myt) + uo2((y — Myt)^ + 

+ V 2 i{x - Mxtf{y - Myt) + ui 2 (x - Mxt)((y - Mytf + 1^) 

+ ^^22((x — Mxt)^{{y — Myt)^ + 1^) + i^/6) 

+ “n(^^/2) 

+ U2i{x — Mxt)t^ + ui2(y — Myt)t^ 

+ «22(2(x — Mxt){y — Myt)t^) 

+ Poi(-<) 

+ Pn(-(a; - Mxt)t) + po2(-2(y - Myt)t) 

+ P2i(-(x - Mxtft - tV3) +Pi2(-2(x - Mxt){y - Myt)t) 

+ P22(~2((x — MxtY + t^/3)(y — Myt)t), 
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and for p by 


pa{x, y, t) = Poo 

+ pio(i - MJ) + poi (y - Myt) 

+ p 2 o((x - Mrtf + i") + pii(x - Mj){y - Myt) +po2{{y - Myt)^ + 1^) 
+ p 2 i{{x - M^tf + t^){y - Myt) + pu{x - M^t){{y - Mytf + 1^) 

+ P22(((x - M,tf + e){{y - Mytf + f2) - 2<V3) 

+ Ulo(~0 

+ u2o(~2(x — Mxt)t) + (y — Myt)t) 

+ U2i(-2(a: - Mxt){y - Myt)t) + «i2(-(y - Mytft - t^S) 

+ U22{-2{x - Mxt){{y - Mytf + ^t'^)t) 

+ Voi{—t) 

+ vii(— (x — Mxt)t) + uo2(— 2(y — Myt)t) 

+ V2i{-{x - Mxtft - t®/3) + ui2(-2(x - Mxt){y - Myt)t) 

+ i;22(-2((x - Mxtf + <V3)(y - Myt)t). 



Appendix D: A Fourth Order Algorithm For Equation (31) 
(Appendix D-U) u solution for the linearized Euler equations 


+k{-pio - MyUoi - MxUio) 

+k^{MyPii + 2MxP2o + MyUo2 + MxMyUii + (1 + M^)u2q + uii/2) 

+k\-{l + 3M2)pi2/3 - 2MxMyP2i - (1 + 3M|)p3o 

- MJuo3 - MxM^Ui2 - (1 + Ml)MyU2l - Mx{Z + MI)u3q 

- MyVi2 - MxV2i) 

+ k\My(l + M2)pi3 + 2Mx{l/Z + MI)p22 
+ (1 + 3iW^)iV/yP3i + 4Afi(l + M^')p4o 
+ MyUQ4 + MxMyU\3 + (1/6 + My 4” M^My')U22 

+ Mx{Z + Ml)MyU3x + (1 + QMl + M*)u4o 
+ (1 + 6M2)ui3/4 + 2MxMyV22 + (1 + 6M2)u3i/4) 

+k\-{l + lOM^ + 5M^)pul5 - 2MxMy(l + M^)p23 

- (1 + 5M^ + 5Mj + 15M2 mJ)p32/5 - 4M*(1 + Ml)MyP4i 

- MxM^ui4 - My(l/2 + + MlMl)u2z 

- M*(l/2 + ZMl + MlMl)u32 - (1 + 6M^ + Mt)MyU4i 

- My{\ + 2MJ)ui 4 - Mx(l/2 + 3 MJ)u23 

- (1/2 + ZMl)MyVz2 - Mxil + 2Ml)v4i) 

+Jfc®((2Mj5 + AMxMl + 2MxM*)p24 

+ (AMx/5 4- AMl/Z + AMxMy + 4M^My)p42 
+ (1/15 + + Mj + MIM*)u24 

+ (2/15 + Ml + Ml + mlMl + MtMl)u42 
+ {2MxMy + 4 MxMI)v24 
+ (2MxMy + 4MlMy)v42) 
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(Appendix D-V) v solution for the linearized Euler equations 


+A:(— poi — MyVoi — MxVio) 

+P(2MyPo2 + MxPii + uii/2 + (1 + My)vo2 + MxMyVii + M^U2o) 

+ 3M^)po3 - 2MxMyPi2 - (1 + 3M2)p2i/3 
— MyU\2 — MxU21 

- My(3 + MI)vo 3 - Mx{l + M^)vi 2 - MlMyV2i - M^uao) 
+k*{4:My{l + My)po4 + Mx(l + ZMy)p\z 

+ 2(1/3 + Ml)MyP22 + Mx{\ + Ml)pii 
+ (1 + &My)uizl^ + 2MxMyU22 + (1 + 6Mj)u3i/4 
+ (1 + 6Mj + M^)uo 4 + MxMy{3 + M^)vi 3 
+ ( 1/6 + + MlMy)v22 + MlMyVzi + MIvaq) 

■irk^{~AMxMy{\ + MJ)pi4 - (1 + bMl + 5MJ + 15M2mJ)p23/5 

- 2Mx{l + Ml)MyP32 - (1 + lOM^ + 5M^)p4i/5 

- My(l + 2MJ)ui 4 - M*(l/2 + 3 MJ)u23 

- (1/2 + 3M^)MyUz2 — Afi(l + 2Afj)«4i 

- Mx(l + 6Mj + Mj)uu - Mj,(l/2 + 3M| + MIM^)v23 

- Mx(l/2 + M^+ MImI)v 32 - M^MyV^i) 

+jfc®((4M„/5 + AMlMy + 4MJ/3 + 4MIM^)p24 

+ (2Mj,/5 + AMlMy + 2M*My)p42 

+ {2MxMy + AMxMy)u24 
+ {2MxMy + AMlMy)u42 

+ (2/15 + + Mj + mlMl + MIM^)v24 

+ (1/15 A- Ml + Ml A- MImI)v42) 
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(Appendix D-P) p solution for the linearized Etiler equations 


+k{-{MyPoi) - MxPio - Uio - Voi) 

+ My)po2 + MxMyPii + (1 + M^)p2o 

+ MyUii + 2MxU20 + 2MyV02 + MxVn) 

+k^{—My{Z + My)po3 — Mi(l + My)pi2 

- (1 + Ml)MyP2x - Mx(Z + M2)p30 

- (1 + 3M2)ui2/3 - 2MxMyU2i - (1 + ZM^)u3o 

- (1 + ZM^)vo3 - 2MxMyVi2 - (1 + 3Ml)v2ilZ) 

+k\(l + ml + Ml)po^ + MxMy{Z + MJ)pi3 

+ (1/3 + + MImI)p22 + Mx{Z + Ml)MyP3\ 

+ (1 + 6M"+M,")p40 
+ My{\ + Ml)ui3 + 2Mx(l/Z + MI)u22 
+ (1 + ZMl')MyU3i + 4Afj;(l + Mj.')v,4q 
+ 4My{l + My )uo4 + Mj(l + 3My )ui3 

+ 2(1/3 + Ml)MyV22 + Mi(l + Mj)u3i) 

+k\-Mx{l + ml + MJ)pi 4 - Mj,(l + ZMl + Ml + MlMl)p2z 

- Mx{l + m2 + ZMl + MlMl)p32 - (1 + ml + Ml)MyP4i 

- (1 + 10M2 + 5Ml)uij5 - 2MxMy(l + MJ)«23 

- (1 + 5Ml + ml + 15MlMl)u32/5 - 4Mx(l + M2 )Mj,U4i 

- AMxMyil + M2)vi4 - (1 + 5M2 + 5M2 + 15M2M2)u23/5 

- 2Mx(l + Ml)MyV32 — (1 + lOMl + 5Mj)u4i/5) 

+fc®((l/5 + Ml + 2Ml + 6MlMl + MJ + MImI)p24 

+ (1/5 + 2Ml + Ml + Ml+ mlMl + MlMl)p42 
+ (2Mx /5 + 4M^M2 + 2MxMI)u24 
+ (4MJ5 + AMUZ + AMxMl + AMImI)u^2 
+ (4Mj,/5 + AMlMy + 4MJ/3 + AMImI)v24 
+ (2My/5 + AMlMy + 2MlMy)v42) 
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(a) Stencils for second, fourth and sixth order interpolation. 
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(b) Square stencil (x,y) expansion coefficients. 
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(c) Minimal stencil (x,y) expansion coefficients. 
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(d) Taylor series (x,y) expansion coefficients. 

Figure 1 . — ^Two dimensional interpolation information. 
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Figure 2. — Plots for the linearized Euler equations with = 0.6 and My = 0, using the fourth order exact propagator algorithm 
along a wall, (a) Pressure at t = 1 5. (b) Pressure at t = 45. (c) Pressure at t = 75. (d) Pressure at t = 1 50. 
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Figure 3. — Plots for the linearized Euler equations with Mjj = 0.5 and My = 0, using the fourth order exact propagator algorithm 
in a duct, (a) Pressure at t = 15. (b) Pressure at t = 45. (c) Pressure at t = 75. (d) Pressure at t = 1 50. 
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